#include "kernels.cuh"

__global__  void multiplyKernel_8_noPrefetch_granularity_1x1_noLoopUnrolling(
	float *result, const float *matrix1, const float *matrix2, const int height, const int width, const int length) { 
		__shared__ float mArray[8][8]; 
		__shared__ float nArray[8][8];

		int bx = blockIdx.x; 
		int by = blockIdx.y; 
		int tx = threadIdx.x; 
		int ty = threadIdx.y; 

		int row = by * 8 + ty;
		int col = bx * 8 + tx;

		float value = 0;
		for (int m = 0; m < (length + 7) / 8; ++m) {
			mArray[ty][tx] = (((m * 8 + tx) < length) && (row < height))? matrix1[row * length + m * 8 + tx] : 0;
			nArray[ty][tx] = (((m * 8 + ty) < length) && (col < width))? matrix2[(m * 8 + ty) * width + col] : 0;
			__syncthreads();
			for (int k = 0; k < 8; ++k) {
				value += mArray[ty][k] * nArray[k][tx]; 
			}
			__syncthreads();
		}
		if (row < height && col < width)
			result[row * width + col] = value;
		
}

__global__  void multiplyKernel_8_noPrefetch_granularity_1x1_2LoopUnrolling(
	float *result, const float *matrix1, const float *matrix2, const int height, const int width, const int length) {
		__shared__ float mArray[8][8]; 
		__shared__ float nArray[8][8];

		int bx = blockIdx.x; 
		int by = blockIdx.y; 
		int tx = threadIdx.x; 
		int ty = threadIdx.y; 

		int row = by * 8 + ty;
		int col = bx * 8 + tx;

		float value = 0;
		for (int m = 0; m < (length + 7) / 8; ++m) {
			mArray[ty][tx] = (((m * 8 + tx) < length) && (row < height))? matrix1[row * length + m * 8 + tx] : 0;
			nArray[ty][tx] = (((m * 8 + ty) < length) && (col < width))? matrix2[(m * 8 + ty) * width + col] : 0;
			__syncthreads();
			for (int k = 0; k < 4; ++k) {
				value += mArray[ty][2 * k] * nArray[2 * k][tx] + 
					mArray[ty][2 * k + 1] * nArray[2 * k + 1][tx]; 
			}
			__syncthreads();
		}
		if (row < height && col < width)
			result[row * width + col] = value;

}

__global__  void multiplyKernel_8_noPrefetch_granularity_1x1_4LoopUnrolling(
	float *result, const float *matrix1, const float *matrix2, const int height, const int width, const int length) { 
		__shared__ float mArray[8][8]; 
		__shared__ float nArray[8][8];

		int bx = blockIdx.x; 
		int by = blockIdx.y; 
		int tx = threadIdx.x; 
		int ty = threadIdx.y; 

		int row = by * 8 + ty;
		int col = bx * 8 + tx;

		float value = 0;
		for (int m = 0; m < (length + 7) / 8; ++m) {
			mArray[ty][tx] = (((m * 8 + tx) < length) && (row < height))? matrix1[row * length + m * 8 + tx] : 0;
			nArray[ty][tx] = (((m * 8 + ty) < length) && (col < width))? matrix2[(m * 8 + ty) * width + col] : 0;
			__syncthreads();
			for (int k = 0; k < 2; ++k) {
				value += mArray[ty][4 * k] * nArray[4 * k][tx] +
					mArray[ty][4 * k + 1] * nArray[4 * k + 1][tx] +
					mArray[ty][4 * k + 2] * nArray[4 * k + 2][tx] +
					mArray[ty][4 * k + 3] * nArray[4 * k + 3][tx]; 
			}
			__syncthreads();
		}
		if (row < height && col < width)
			result[row * width + col] = value;
		
}

__global__  void multiplyKernel_8_noPrefetch_granularity_1x1_fullLoopUnrolling(
	float *result, const float *matrix1, const float *matrix2, const int height, const int width, const int length) { 
		__shared__ float mArray[8][8]; 
		__shared__ float nArray[8][8];

		int bx = blockIdx.x; 
		int by = blockIdx.y; 
		int tx = threadIdx.x; 
		int ty = threadIdx.y; 

		int row = by * 8 + ty;
		int col = bx * 8 + tx;

		float value = 0;
		for (int m = 0; m < (length + 7) / 8; ++m) {
			mArray[ty][tx] = (((m * 8 + tx) < length) && (row < height))? matrix1[row * length + m * 8 + tx] : 0;
			nArray[ty][tx] = (((m * 8 + ty) < length) && (col < width))? matrix2[(m * 8 + ty) * width + col] : 0;
			__syncthreads();
			value += mArray[ty][0] * nArray[0][tx] +
				mArray[ty][1] * nArray[1][tx] +
				mArray[ty][2] * nArray[2][tx] +
				mArray[ty][3] * nArray[3][tx] +
				mArray[ty][4] * nArray[4][tx] +
				mArray[ty][5] * nArray[5][tx] +
				mArray[ty][6] * nArray[6][tx] +
				mArray[ty][7] * nArray[7][tx];
			__syncthreads();
		}
		if (row < height && col < width)
			result[row * width + col] = value;
}

__global__  void multiplyKernel_8_noPrefetch_granularity_1x2_noLoopUnrolling(
	float *result, const float *matrix1, const float *matrix2, const int height, const int width, const int length) { 
		__shared__ float mArray[8][8]; 
		__shared__ float nArray0[8][8];
		__shared__ float nArray1[8][8];
		
		int bx = blockIdx.x; 
		int by = blockIdx.y; 
		int tx = threadIdx.x; 
		int ty = threadIdx.y; 

		int row = by * 8 + ty;
		int col0 = 2 * bx * 8 + tx;
		int col1 = (2 * bx + 1) * 8 + tx;

		float value0 = 0;
		float value1 = 0;
		for (int m = 0; m < (length + 7) / 8; ++m) {
			mArray[ty][tx] = (((m * 8 + tx) < length) && (row < height))? matrix1[row * length + m * 8 + tx] : 0;
			nArray0[ty][tx] = (((m * 8 + ty) < length) && (col0 < width))? matrix2[(m * 8 + ty) * width + col0] : 0;
			nArray1[ty][tx] = (((m * 8 + ty) < length) && (col1 < width))? matrix2[(m * 8 + ty) * width + col1] : 0;
			__syncthreads();
			for (int k = 0; k < 8; ++k) {
				value0 += mArray[ty][k] * nArray0[k][tx];
				value1 += mArray[ty][k] * nArray1[k][tx];
			}
			__syncthreads();
		}
		if (row < height) {
			if (col0 < width)
				result[row * width + col0] = value0;
			if (col1 < width)
				result[row * width + col1] = value1;
		}
}

__global__  void multiplyKernel_8_noPrefetch_granularity_1x2_2LoopUnrolling(
	float *result, const float *matrix1, const float *matrix2, const int height, const int width, const int length) { 
		__shared__ float mArray[8][8]; 
		__shared__ float nArray0[8][8];
		__shared__ float nArray1[8][8];
		
		int bx = blockIdx.x; 
		int by = blockIdx.y; 
		int tx = threadIdx.x; 
		int ty = threadIdx.y; 

		int row = by * 8 + ty;
		int col0 = 2 * bx * 8 + tx;
		int col1 = (2 * bx + 1) * 8 + tx;

		float value0 = 0;
		float value1 = 0;
		for (int m = 0; m < (length + 7) / 8; ++m) {
			mArray[ty][tx] = (((m * 8 + tx) < length) && (row < height))? matrix1[row * length + m * 8 + tx] : 0;
			nArray0[ty][tx] = (((m * 8 + ty) < length) && (col0 < width))? matrix2[(m * 8 + ty) * width + col0] : 0;
			nArray1[ty][tx] = (((m * 8 + ty) < length) && (col1 < width))? matrix2[(m * 8 + ty) * width + col1] : 0;
			__syncthreads();
			for (int k = 0; k < 4; ++k) {
				value0 += mArray[ty][2 * k] * nArray0[2 * k][tx] + 
					mArray[ty][2 * k + 1] * nArray0[2 * k + 1][tx]; 
				value1 += mArray[ty][2 * k] * nArray1[2 * k][tx] + 
					mArray[ty][2 * k + 1] * nArray1[2 * k + 1][tx]; 
			}
			__syncthreads();
		}
		if (row < height) {
			if (col0 < width)
				result[row * width + col0] = value0;
			if (col1 < width)
				result[row * width + col1] = value1;
		}
}

__global__  void multiplyKernel_8_noPrefetch_granularity_1x2_4LoopUnrolling(
	float *result, const float *matrix1, const float *matrix2, const int height, const int width, const int length) { 
		__shared__ float mArray[8][8]; 
		__shared__ float nArray0[8][8];
		__shared__ float nArray1[8][8];
		
		int bx = blockIdx.x; 
		int by = blockIdx.y; 
		int tx = threadIdx.x; 
		int ty = threadIdx.y; 

		int row = by * 8 + ty;
		int col0 = 2 * bx * 8 + tx;
		int col1 = (2 * bx + 1) * 8 + tx;

		float value0 = 0;
		float value1 = 0;
		for (int m = 0; m < (length + 7) / 8; ++m) {
			mArray[ty][tx] = (((m * 8 + tx) < length) && (row < height))? matrix1[row * length + m * 8 + tx] : 0;
			nArray0[ty][tx] = (((m * 8 + ty) < length) && (col0 < width))? matrix2[(m * 8 + ty) * width + col0] : 0;
			nArray1[ty][tx] = (((m * 8 + ty) < length) && (col1 < width))? matrix2[(m * 8 + ty) * width + col1] : 0;
			__syncthreads();
			for (int k = 0; k < 2; ++k) {
				value0 += mArray[ty][4 * k] * nArray0[4 * k][tx] +
					mArray[ty][4 * k + 1] * nArray0[4 * k + 1][tx] +
					mArray[ty][4 * k + 2] * nArray0[4 * k + 2][tx] +
					mArray[ty][4 * k + 3] * nArray0[4 * k + 3][tx]; 
				value1 += mArray[ty][4 * k] * nArray1[4 * k][tx] +
					mArray[ty][4 * k + 1] * nArray1[4 * k + 1][tx] +
					mArray[ty][4 * k + 2] * nArray1[4 * k + 2][tx] +
					mArray[ty][4 * k + 3] * nArray1[4 * k + 3][tx];		
			}
			__syncthreads();
		}
		if (row < height) {
			if (col0 < width)
				result[row * width + col0] = value0;
			if (col1 < width)
				result[row * width + col1] = value1;
		}
}

__global__  void multiplyKernel_8_noPrefetch_granularity_1x2_fullLoopUnrolling(
	float *result, const float *matrix1, const float *matrix2, const int height, const int width, const int length) { 
		__shared__ float mArray[8][8]; 
		__shared__ float nArray0[8][8];
		__shared__ float nArray1[8][8];
		
		int bx = blockIdx.x; 
		int by = blockIdx.y; 
		int tx = threadIdx.x; 
		int ty = threadIdx.y; 

		int row = by * 8 + ty;
		int col0 = 2 * bx * 8 + tx;
		int col1 = (2 * bx + 1) * 8 + tx;

		float value0 = 0;
		float value1 = 0;
		for (int m = 0; m < (length + 7) / 8; ++m) {
			mArray[ty][tx] = (((m * 8 + tx) < length) && (row < height))? matrix1[row * length + m * 8 + tx] : 0;
			nArray0[ty][tx] = (((m * 8 + ty) < length) && (col0 < width))? matrix2[(m * 8 + ty) * width + col0] : 0;
			nArray1[ty][tx] = (((m * 8 + ty) < length) && (col1 < width))? matrix2[(m * 8 + ty) * width + col1] : 0;
			__syncthreads();
			value0 += mArray[ty][0] * nArray0[0][tx] +
				mArray[ty][1] * nArray0[1][tx] +
				mArray[ty][2] * nArray0[2][tx] +
				mArray[ty][3] * nArray0[3][tx] +
				mArray[ty][4] * nArray0[4][tx] +
				mArray[ty][5] * nArray0[5][tx] +
				mArray[ty][6] * nArray0[6][tx] +
				mArray[ty][7] * nArray0[7][tx];
			value1 += mArray[ty][0] * nArray1[0][tx] +
				mArray[ty][1] * nArray1[1][tx] +
				mArray[ty][2] * nArray1[2][tx] +
				mArray[ty][3] * nArray1[3][tx] +
				mArray[ty][4] * nArray1[4][tx] +
				mArray[ty][5] * nArray1[5][tx] +
				mArray[ty][6] * nArray1[6][tx] +
				mArray[ty][7] * nArray1[7][tx];
			__syncthreads();
		}
		if (row < height) {
			if (col0 < width)
				result[row * width + col0] = value0;
			if (col1 < width)
				result[row * width + col1] = value1;
		}
}

__global__  void multiplyKernel_8_noPrefetch_granularity_1x4_noLoopUnrolling(
	float *result, const float *matrix1, const float *matrix2, const int height, const int width, const int length) { 
		__shared__ float mArray[8][8]; 
		__shared__ float nArray0[8][8];
		__shared__ float nArray1[8][8];
		__shared__ float nArray2[8][8];
		__shared__ float nArray3[8][8];
		
		int bx = blockIdx.x; 
		int by = blockIdx.y; 
		int tx = threadIdx.x; 
		int ty = threadIdx.y; 

		int row = by * 8 + ty;
		int col0 = 4 * bx * 8 + tx;
		int col1 = (4 * bx + 1) * 8 + tx;
		int col2 = (4 * bx + 2) * 8 + tx;
		int col3 = (4 * bx + 3) * 8 + tx;

		float value0 = 0;
		float value1 = 0;
		float value2 = 0;
		float value3 = 0;
		for (int m = 0; m < (length + 7) / 8; ++m) {
			mArray[ty][tx] = (((m * 8 + tx) < length) && (row < height))? matrix1[row * length + m * 8 + tx] : 0;
			nArray0[ty][tx] = (((m * 8 + ty) < length) && (col0 < width))? matrix2[(m * 8 + ty) * width + col0] : 0;
			nArray1[ty][tx] = (((m * 8 + ty) < length) && (col1 < width))? matrix2[(m * 8 + ty) * width + col1] : 0;
			nArray2[ty][tx] = (((m * 8 + ty) < length) && (col2 < width))? matrix2[(m * 8 + ty) * width + col2] : 0;
			nArray3[ty][tx] = (((m * 8 + ty) < length) && (col3 < width))? matrix2[(m * 8 + ty) * width + col3] : 0;
			__syncthreads();
			for (int k = 0; k < 8; ++k) {
				value0 += mArray[ty][k] * nArray0[k][tx];
				value1 += mArray[ty][k] * nArray1[k][tx];
				value2 += mArray[ty][k] * nArray2[k][tx];
				value3 += mArray[ty][k] * nArray3[k][tx];
			}
			__syncthreads();
		}
		if (row < height) {
			if (col0 < width)
				result[row * width + col0] = value0;
			if (col1 < width)
				result[row * width + col1] = value1;
			if (col2 < width)
				result[row * width + col2] = value2;
			if (col3 < width)
				result[row * width + col3] = value3;
		}
}

__global__  void multiplyKernel_8_noPrefetch_granularity_1x4_2LoopUnrolling(
	float *result, const float *matrix1, const float *matrix2, const int height, const int width, const int length) { 
		__shared__ float mArray[8][8]; 
		__shared__ float nArray0[8][8];
		__shared__ float nArray1[8][8];
		__shared__ float nArray2[8][8];
		__shared__ float nArray3[8][8];
		
		int bx = blockIdx.x; 
		int by = blockIdx.y; 
		int tx = threadIdx.x; 
		int ty = threadIdx.y; 

		int row = by * 8 + ty;
		int col0 = 4 * bx * 8 + tx;
		int col1 = (4 * bx + 1) * 8 + tx;
		int col2 = (4 * bx + 2) * 8 + tx;
		int col3 = (4 * bx + 3) * 8 + tx;

		float value0 = 0;
		float value1 = 0;
		float value2 = 0;
		float value3 = 0;
		for (int m = 0; m < (length + 7) / 8; ++m) {
			mArray[ty][tx] = (((m * 8 + tx) < length) && (row < height))? matrix1[row * length + m * 8 + tx] : 0;
			nArray0[ty][tx] = (((m * 8 + ty) < length) && (col0 < width))? matrix2[(m * 8 + ty) * width + col0] : 0;
			nArray1[ty][tx] = (((m * 8 + ty) < length) && (col1 < width))? matrix2[(m * 8 + ty) * width + col1] : 0;
			nArray2[ty][tx] = (((m * 8 + ty) < length) && (col2 < width))? matrix2[(m * 8 + ty) * width + col2] : 0;
			nArray3[ty][tx] = (((m * 8 + ty) < length) && (col3 < width))? matrix2[(m * 8 + ty) * width + col3] : 0;
			__syncthreads();
			for (int k = 0; k < 4; ++k) {
				value0 += mArray[ty][2 * k] * nArray0[2 * k][tx] + 
					mArray[ty][2 * k + 1] * nArray0[2 * k + 1][tx];
				value1 += mArray[ty][2 * k] * nArray1[2 * k][tx] + 
					mArray[ty][2 * k + 1] * nArray1[2 * k + 1][tx];
				value2 += mArray[ty][2 * k] * nArray2[2 * k][tx] + 
					mArray[ty][2 * k + 1] * nArray2[2 * k + 1][tx];
				value3 += mArray[ty][2 * k] * nArray3[2 * k][tx] + 
					mArray[ty][2 * k + 1] * nArray3[2 * k + 1][tx];
			}
			__syncthreads();
		}
		if (row < height) {
			if (col0 < width)
				result[row * width + col0] = value0;
			if (col1 < width)
				result[row * width + col1] = value1;
			if (col2 < width)
				result[row * width + col2] = value2;
			if (col3 < width)
				result[row * width + col3] = value3;
		}
}

__global__  void multiplyKernel_8_noPrefetch_granularity_1x4_4LoopUnrolling(
	float *result, const float *matrix1, const float *matrix2, const int height, const int width, const int length) { 
		__shared__ float mArray[8][8]; 
		__shared__ float nArray0[8][8];
		__shared__ float nArray1[8][8];
		__shared__ float nArray2[8][8];
		__shared__ float nArray3[8][8];
		
		int bx = blockIdx.x; 
		int by = blockIdx.y; 
		int tx = threadIdx.x; 
		int ty = threadIdx.y; 

		int row = by * 8 + ty;
		int col0 = 4 * bx * 8 + tx;
		int col1 = (4 * bx + 1) * 8 + tx;
		int col2 = (4 * bx + 2) * 8 + tx;
		int col3 = (4 * bx + 3) * 8 + tx;

		float value0 = 0;
		float value1 = 0;
		float value2 = 0;
		float value3 = 0;
		for (int m = 0; m < (length + 7) / 8; ++m) {
			mArray[ty][tx] = (((m * 8 + tx) < length) && (row < height))? matrix1[row * length + m * 8 + tx] : 0;
			nArray0[ty][tx] = (((m * 8 + ty) < length) && (col0 < width))? matrix2[(m * 8 + ty) * width + col0] : 0;
			nArray1[ty][tx] = (((m * 8 + ty) < length) && (col1 < width))? matrix2[(m * 8 + ty) * width + col1] : 0;
			nArray2[ty][tx] = (((m * 8 + ty) < length) && (col2 < width))? matrix2[(m * 8 + ty) * width + col2] : 0;
			nArray3[ty][tx] = (((m * 8 + ty) < length) && (col3 < width))? matrix2[(m * 8 + ty) * width + col3] : 0;
			__syncthreads();
			for (int k = 0; k < 2; ++k) {
				value0 += mArray[ty][4 * k] * nArray0[4 * k][tx] + 
					mArray[ty][4 * k + 1] * nArray0[4 * k + 1][tx] + 
					mArray[ty][4 * k + 2] * nArray0[4 * k + 2][tx] + 
					mArray[ty][4 * k + 3] * nArray0[4 * k + 3][tx];
				value1 += mArray[ty][4 * k] * nArray1[4 * k][tx] + 
					mArray[ty][4 * k + 1] * nArray1[4 * k + 1][tx] + 
					mArray[ty][4 * k + 2] * nArray1[4 * k + 2][tx] + 
					mArray[ty][4 * k + 3] * nArray1[4 * k + 3][tx];
				value2 += mArray[ty][4 * k] * nArray2[4 * k][tx] + 
					mArray[ty][4 * k + 1] * nArray2[4 * k + 1][tx] + 
					mArray[ty][4 * k + 2] * nArray2[4 * k + 2][tx] + 
					mArray[ty][4 * k + 3] * nArray2[4 * k + 3][tx];
				value3 += mArray[ty][4 * k] * nArray3[4 * k][tx] + 
					mArray[ty][4 * k + 1] * nArray3[4 * k + 1][tx] + 
					mArray[ty][4 * k + 2] * nArray3[4 * k + 2][tx] + 
					mArray[ty][4 * k + 3] * nArray3[4 * k + 3][tx];
			}
			__syncthreads();
		}
		if (row < height) {
			if (col0 < width)
				result[row * width + col0] = value0;
			if (col1 < width)
				result[row * width + col1] = value1;
			if (col2 < width)
				result[row * width + col2] = value2;
			if (col3 < width)
				result[row * width + col3] = value3;
		}
}

__global__  void multiplyKernel_8_noPrefetch_granularity_1x4_fullLoopUnrolling(
	float *result, const float *matrix1, const float *matrix2, const int height, const int width, const int length) { 
		__shared__ float mArray[8][8]; 
		__shared__ float nArray0[8][8];
		__shared__ float nArray1[8][8];
		__shared__ float nArray2[8][8];
		__shared__ float nArray3[8][8];
		
		int bx = blockIdx.x; 
		int by = blockIdx.y; 
		int tx = threadIdx.x; 
		int ty = threadIdx.y; 

		int row = by * 8 + ty;
		int col0 = 4 * bx * 8 + tx;
		int col1 = (4 * bx + 1) * 8 + tx;
		int col2 = (4 * bx + 2) * 8 + tx;
		int col3 = (4 * bx + 3) * 8 + tx;

		float value0 = 0;
		float value1 = 0;
		float value2 = 0;
		float value3 = 0;
		for (int m = 0; m < (length + 7) / 8; ++m) {
			mArray[ty][tx] = (((m * 8 + tx) < length) && (row < height))? matrix1[row * length + m * 8 + tx] : 0;
			nArray0[ty][tx] = (((m * 8 + ty) < length) && (col0 < width))? matrix2[(m * 8 + ty) * width + col0] : 0;
			nArray1[ty][tx] = (((m * 8 + ty) < length) && (col1 < width))? matrix2[(m * 8 + ty) * width + col1] : 0;
			nArray2[ty][tx] = (((m * 8 + ty) < length) && (col2 < width))? matrix2[(m * 8 + ty) * width + col2] : 0;
			nArray3[ty][tx] = (((m * 8 + ty) < length) && (col3 < width))? matrix2[(m * 8 + ty) * width + col3] : 0;
			__syncthreads();

			value0 += mArray[ty][0] * nArray0[0][tx] +
				mArray[ty][1] * nArray0[1][tx] +
				mArray[ty][2] * nArray0[2][tx] +
				mArray[ty][3] * nArray0[3][tx] +
				mArray[ty][4] * nArray0[4][tx] +
				mArray[ty][5] * nArray0[5][tx] +
				mArray[ty][6] * nArray0[6][tx] +
				mArray[ty][7] * nArray0[7][tx];
			value1 += mArray[ty][0] * nArray1[0][tx] +
				mArray[ty][1] * nArray1[1][tx] +
				mArray[ty][2] * nArray1[2][tx] +
				mArray[ty][3] * nArray1[3][tx] +
				mArray[ty][4] * nArray1[4][tx] +
				mArray[ty][5] * nArray1[5][tx] +
				mArray[ty][6] * nArray1[6][tx] +
				mArray[ty][7] * nArray1[7][tx];
			value2 += mArray[ty][0] * nArray2[0][tx] +
				mArray[ty][1] * nArray2[1][tx] +
				mArray[ty][2] * nArray2[2][tx] +
				mArray[ty][3] * nArray2[3][tx] +
				mArray[ty][4] * nArray2[4][tx] +
				mArray[ty][5] * nArray2[5][tx] +
				mArray[ty][6] * nArray2[6][tx] +
				mArray[ty][7] * nArray2[7][tx];
			value3 += mArray[ty][0] * nArray3[0][tx] +
				mArray[ty][1] * nArray3[1][tx] +
				mArray[ty][2] * nArray3[2][tx] +
				mArray[ty][3] * nArray3[3][tx] +
				mArray[ty][4] * nArray3[4][tx] +
				mArray[ty][5] * nArray3[5][tx] +
				mArray[ty][6] * nArray3[6][tx] +
				mArray[ty][7] * nArray3[7][tx];
			__syncthreads();
		}
		if (row < height) {
			if (col0 < width)
				result[row * width + col0] = value0;
			if (col1 < width)
				result[row * width + col1] = value1;
			if (col2 < width)
				result[row * width + col2] = value2;
			if (col3 < width)
				result[row * width + col3] = value3;
		}
}

__global__  void multiplyKernel_8_prefetch_granularity_1x1_noLoopUnrolling(
	float *result, const float *matrix1, const float *matrix2, const int height, const int width, const int length) { 
		__shared__ float mArray[8][8]; 
		__shared__ float nArray[8][8];

		int bx = blockIdx.x; 
		int by = blockIdx.y; 
		int tx = threadIdx.x; 
		int ty = threadIdx.y; 

		int row = by * 8 + ty;
		int col = bx * 8 + tx;

		float value = 0;
		float mElement = ((tx < length) && (row < height))? matrix1[row * length + tx] : 0;
		float nElement = ((ty < length) && (col < width))? matrix2[ty * width + col] : 0;
		for (int m = 1; m <= (length + 7) / 8; ++m) {
			mArray[ty][tx] = mElement;
			nArray[ty][tx] = nElement;
			__syncthreads();
			mElement = (((m * 8 + tx) < length) && (row < height))? matrix1[row * length + m * 8 + tx] : 0;
			nElement = (((m * 8 + ty) < length) && (col < width))? matrix2[(m * 8 + ty) * width + col] : 0;
			for (int k = 0; k < 8; ++k) {
				value += mArray[ty][k] * nArray[k][tx]; 
			}
			__syncthreads();
		}
		if (row < height && col < width)
			result[row * width + col] = value;
}

__global__  void multiplyKernel_8_prefetch_granularity_1x1_2LoopUnrolling(
	float *result, const float *matrix1, const float *matrix2, const int height, const int width, const int length) { 
		__shared__ float mArray[8][8]; 
		__shared__ float nArray[8][8];

		int bx = blockIdx.x; 
		int by = blockIdx.y; 
		int tx = threadIdx.x; 
		int ty = threadIdx.y; 

		int row = by * 8 + ty;
		int col = bx * 8 + tx;

		float value = 0;
		float mElement = ((tx < length) && (row < height))? matrix1[row * length + tx] : 0;
		float nElement = ((ty < length) && (col < width))? matrix2[ty * width + col] : 0;
		for (int m = 1; m <= (length + 7) / 8; ++m) {
			mArray[ty][tx] = mElement;
			nArray[ty][tx] = nElement;
			__syncthreads();
			mElement = (((m * 8 + tx) < length) && (row < height))? matrix1[row * length + m * 8 + tx] : 0;
			nElement = (((m * 8 + ty) < length) && (col < width))? matrix2[(m * 8 + ty) * width + col] : 0;
			for (int k = 0; k < 4; ++k) {
				value += mArray[ty][2 * k] * nArray[2 * k][tx] + 
					mArray[ty][2 * k + 1] * nArray[2 * k + 1][tx]; 
			}
			__syncthreads();
		}
		if (row < height && col < width)
			result[row * width + col] = value;
}

__global__  void multiplyKernel_8_prefetch_granularity_1x1_4LoopUnrolling(
	float *result, const float *matrix1, const float *matrix2, const int height, const int width, const int length) { 
		__shared__ float mArray[8][8]; 
		__shared__ float nArray[8][8];

		int bx = blockIdx.x; 
		int by = blockIdx.y; 
		int tx = threadIdx.x; 
		int ty = threadIdx.y; 

		int row = by * 8 + ty;
		int col = bx * 8 + tx;

		float value = 0;
		float mElement = ((tx < length) && (row < height))? matrix1[row * length + tx] : 0;
		float nElement = ((ty < length) && (col < width))? matrix2[ty * width + col] : 0;
		for (int m = 1; m <= (length + 7) / 8; ++m) {
			mArray[ty][tx] = mElement;
			nArray[ty][tx] = nElement;
			__syncthreads();
			mElement = (((m * 8 + tx) < length) && (row < height))? matrix1[row * length + m * 8 + tx] : 0;
			nElement = (((m * 8 + ty) < length) && (col < width))? matrix2[(m * 8 + ty) * width + col] : 0;
			for (int k = 0; k < 2; ++k) {
				value += mArray[ty][4 * k] * nArray[4 * k][tx] +
					mArray[ty][4 * k + 1] * nArray[4 * k + 1][tx] +
					mArray[ty][4 * k + 2] * nArray[4 * k + 2][tx] +
					mArray[ty][4 * k + 3] * nArray[4 * k + 3][tx]; 
			}
			__syncthreads();
		}
		if (row < height && col < width)
			result[row * width + col] = value;
}

__global__  void multiplyKernel_8_prefetch_granularity_1x1_fullLoopUnrolling(
	float *result, const float *matrix1, const float *matrix2, const int height, const int width, const int length) { 
		__shared__ float mArray[8][8]; 
		__shared__ float nArray[8][8];

		int bx = blockIdx.x; 
		int by = blockIdx.y; 
		int tx = threadIdx.x; 
		int ty = threadIdx.y; 

		int row = by * 8 + ty;
		int col = bx * 8 + tx;

		float value = 0;
		float mElement = ((tx < length) && (row < height))? matrix1[row * length + tx] : 0;
		float nElement = ((ty < length) && (col < width))? matrix2[ty * width + col] : 0;
		for (int m = 1; m <= (length + 7) / 8; ++m) {
			mArray[ty][tx] = mElement;
			nArray[ty][tx] = nElement;
			__syncthreads();
			mElement = (((m * 8 + tx) < length) && (row < height))? matrix1[row * length + m * 8 + tx] : 0;
			nElement = (((m * 8 + ty) < length) && (col < width))? matrix2[(m * 8 + ty) * width + col] : 0;
			value += mArray[ty][0] * nArray[0][tx] + 
				mArray[ty][1] * nArray[1][tx] + 
				mArray[ty][2] * nArray[2][tx] + 
				mArray[ty][3] * nArray[3][tx] + 
				mArray[ty][4] * nArray[4][tx] + 
				mArray[ty][5] * nArray[5][tx] + 
				mArray[ty][6] * nArray[6][tx] + 
				mArray[ty][7] * nArray[7][tx];
			__syncthreads();
		}
		if (row < height && col < width)
			result[row * width + col] = value;
}

__global__  void multiplyKernel_8_prefetch_granularity_1x2_noLoopUnrolling(
	float *result, const float *matrix1, const float *matrix2, const int height, const int width, const int length) { 
		__shared__ float mArray[8][8]; 
		__shared__ float nArray0[8][8];
		__shared__ float nArray1[8][8];
		
		int bx = blockIdx.x; 
		int by = blockIdx.y; 
		int tx = threadIdx.x; 
		int ty = threadIdx.y; 

		int row = by * 8 + ty;
		int col0 = 2 * bx * 8 + tx;
		int col1 = (2 * bx + 1) * 8 + tx;

		float value0 = 0;
		float value1 = 0;
		float mElement = ((tx < length) && (row < height))? matrix1[row * length + tx] : 0;
		float nElement0 = ((ty < length) && (col0 < width))? matrix2[ty * width + col0] : 0;
		float nElement1 = ((ty < length) && (col1 < width))? matrix2[ty * width + col1] : 0;
		for (int m = 1; m <= (length + 7) / 8; ++m) {
			mArray[ty][tx] = mElement;
			nArray0[ty][tx] = nElement0;
			nArray1[ty][tx] = nElement1;
			__syncthreads();

			mElement = (((m * 8 + tx) < length) && (row < height))? matrix1[row * length + m * 8 + tx] : 0;
			nElement0 = (((m * 8 + ty) < length) && (col0 < width))? matrix2[(m * 8 + ty) * width + col0] : 0;
			nElement1 = (((m * 8 + ty) < length) && (col1 < width))? matrix2[(m * 8 + ty) * width + col1] : 0;
			
			for (int k = 0; k < 8; ++k) {
				value0 += mArray[ty][k] * nArray0[k][tx];
				value1 += mArray[ty][k] * nArray1[k][tx];
			}
			__syncthreads();
		}
		if (row < height) {
			if (col0 < width)
				result[row * width + col0] = value0;
			if (col1 < width)
				result[row * width + col1] = value1;
		}
}

__global__  void multiplyKernel_8_prefetch_granularity_1x2_2LoopUnrolling(
	float *result, const float *matrix1, const float *matrix2, const int height, const int width, const int length) { 
		__shared__ float mArray[8][8]; 
		__shared__ float nArray0[8][8];
		__shared__ float nArray1[8][8];
		
		int bx = blockIdx.x; 
		int by = blockIdx.y; 
		int tx = threadIdx.x; 
		int ty = threadIdx.y; 

		int row = by * 8 + ty;
		int col0 = 2 * bx * 8 + tx;
		int col1 = (2 * bx + 1) * 8 + tx;

		float value0 = 0;
		float value1 = 0;

		float mElement = ((tx < length) && (row < height))? matrix1[row * length + tx] : 0;
		float nElement0 = ((ty < length) && (col0 < width))? matrix2[ty * width + col0] : 0;
		float nElement1 = ((ty < length) && (col1 < width))? matrix2[ty * width + col1] : 0;

		for (int m = 1; m <= (length + 7) / 8; ++m) {
			mArray[ty][tx] = mElement;
			nArray0[ty][tx] = nElement0;
			nArray1[ty][tx] = nElement1;
			__syncthreads();

			mElement = (((m * 8 + tx) < length) && (row < height))? matrix1[row * length + m * 8 + tx] : 0;
			nElement0 = (((m * 8 + ty) < length) && (col0 < width))? matrix2[(m * 8 + ty) * width + col0] : 0;
			nElement1 = (((m * 8 + ty) < length) && (col1 < width))? matrix2[(m * 8 + ty) * width + col1] : 0;
			
			for (int k = 0; k < 4; ++k) {
				value0 += mArray[ty][2 * k] * nArray0[2 * k][tx] + 
					mArray[ty][2 * k + 1] * nArray0[2 * k + 1][tx]; 
				value1 += mArray[ty][2 * k] * nArray1[2 * k][tx] + 
					mArray[ty][2 * k + 1] * nArray1[2 * k + 1][tx]; 
			}
			__syncthreads();
		}
		if (row < height) {
			if (col0 < width)
				result[row * width + col0] = value0;
			if (col1 < width)
				result[row * width + col1] = value1;
		}
}

__global__  void multiplyKernel_8_prefetch_granularity_1x2_4LoopUnrolling(
	float *result, const float *matrix1, const float *matrix2, const int height, const int width, const int length) { 
		__shared__ float mArray[8][8]; 
		__shared__ float nArray0[8][8];
		__shared__ float nArray1[8][8];
		
		int bx = blockIdx.x; 
		int by = blockIdx.y; 
		int tx = threadIdx.x; 
		int ty = threadIdx.y; 

		int row = by * 8 + ty;
		int col0 = 2 * bx * 8 + tx;
		int col1 = (2 * bx + 1) * 8 + tx;

		float value0 = 0;
		float value1 = 0;

		float mElement = ((tx < length) && (row < height))? matrix1[row * length + tx] : 0;
		float nElement0 = ((ty < length) && (col0 < width))? matrix2[ty * width + col0] : 0;
		float nElement1 = ((ty < length) && (col1 < width))? matrix2[ty * width + col1] : 0;

		for (int m = 1; m <= (length + 7) / 8; ++m) {
			mArray[ty][tx] = mElement;
			nArray0[ty][tx] = nElement0;
			nArray1[ty][tx] = nElement1;
			__syncthreads();

			mElement = (((m * 8 + tx) < length) && (row < height))? matrix1[row * length + m * 8 + tx] : 0;
			nElement0 = (((m * 8 + ty) < length) && (col0 < width))? matrix2[(m * 8 + ty) * width + col0] : 0;
			nElement1 = (((m * 8 + ty) < length) && (col1 < width))? matrix2[(m * 8 + ty) * width + col1] : 0;
			
			for (int k = 0; k < 2; ++k) {
				value0 += mArray[ty][4 * k] * nArray0[4 * k][tx] +
					mArray[ty][4 * k + 1] * nArray0[4 * k + 1][tx] +
					mArray[ty][4 * k + 2] * nArray0[4 * k + 2][tx] +
					mArray[ty][4 * k + 3] * nArray0[4 * k + 3][tx]; 
				value1 += mArray[ty][4 * k] * nArray1[4 * k][tx] +
					mArray[ty][4 * k + 1] * nArray1[4 * k + 1][tx] +
					mArray[ty][4 * k + 2] * nArray1[4 * k + 2][tx] +
					mArray[ty][4 * k + 3] * nArray1[4 * k + 3][tx];		
			}
			__syncthreads();
		}
		if (row < height) {
			if (col0 < width)
				result[row * width + col0] = value0;
			if (col1 < width)
				result[row * width + col1] = value1;
		}
}

__global__  void multiplyKernel_8_prefetch_granularity_1x2_fullLoopUnrolling(
	float *result, const float *matrix1, const float *matrix2, const int height, const int width, const int length) { 
		__shared__ float mArray[8][8]; 
		__shared__ float nArray0[8][8];
		__shared__ float nArray1[8][8];
		
		int bx = blockIdx.x; 
		int by = blockIdx.y; 
		int tx = threadIdx.x; 
		int ty = threadIdx.y; 

		int row = by * 8 + ty;
		int col0 = 2 * bx * 8 + tx;
		int col1 = (2 * bx + 1) * 8 + tx;

		float value0 = 0;
		float value1 = 0;

		float mElement = ((tx < length) && (row < height))? matrix1[row * length + tx] : 0;
		float nElement0 = ((ty < length) && (col0 < width))? matrix2[ty * width + col0] : 0;
		float nElement1 = ((ty < length) && (col1 < width))? matrix2[ty * width + col1] : 0;

		for (int m = 1; m <= (length + 7) / 8; ++m) {
			mArray[ty][tx] = mElement;
			nArray0[ty][tx] = nElement0;
			nArray1[ty][tx] = nElement1;
			__syncthreads();

			mElement = (((m * 8 + tx) < length) && (row < height))? matrix1[row * length + m * 8 + tx] : 0;
			nElement0 = (((m * 8 + ty) < length) && (col0 < width))? matrix2[(m * 8 + ty) * width + col0] : 0;
			nElement1 = (((m * 8 + ty) < length) && (col1 < width))? matrix2[(m * 8 + ty) * width + col1] : 0;
			
			value0 += mArray[ty][0] * nArray0[0][tx] +
					mArray[ty][1] * nArray0[1][tx] +
					mArray[ty][2] * nArray0[2][tx] +
					mArray[ty][3] * nArray0[3][tx] +
					mArray[ty][4] * nArray0[4][tx] +
					mArray[ty][5] * nArray0[5][tx] +
					mArray[ty][6] * nArray0[6][tx] +
					mArray[ty][7] * nArray0[7][tx];
			value1 += mArray[ty][0] * nArray1[0][tx] +
					mArray[ty][1] * nArray1[1][tx] +
					mArray[ty][2] * nArray1[2][tx] +
					mArray[ty][3] * nArray1[3][tx] +
					mArray[ty][4] * nArray1[4][tx] +
					mArray[ty][5] * nArray1[5][tx] +
					mArray[ty][6] * nArray1[6][tx] +
					mArray[ty][7] * nArray1[7][tx];
							
			__syncthreads();
		}
		if (row < height) {
			if (col0 < width)
				result[row * width + col0] = value0;
			if (col1 < width)
				result[row * width + col1] = value1;
		}
}

__global__  void multiplyKernel_8_prefetch_granularity_1x4_noLoopUnrolling(
	float *result, const float *matrix1, const float *matrix2, const int height, const int width, const int length) { 
		__shared__ float mArray[8][8]; 
		__shared__ float nArray0[8][8];
		__shared__ float nArray1[8][8];
		__shared__ float nArray2[8][8];
		__shared__ float nArray3[8][8];
		
		int bx = blockIdx.x; 
		int by = blockIdx.y; 
		int tx = threadIdx.x; 
		int ty = threadIdx.y; 

		int row = by * 8 + ty;
		int col0 = 4 * bx * 8 + tx;
		int col1 = (4 * bx + 1) * 8 + tx;
		int col2 = (4 * bx + 2) * 8 + tx;
		int col3 = (4 * bx + 3) * 8 + tx;

		float value0 = 0;
		float value1 = 0;
		float value2 = 0;
		float value3 = 0;

		float mElement = ((tx < length) && (row < height))? matrix1[row * length + tx] : 0;
		float nElement0 = ((ty < length) && (col0 < width))? matrix2[ty * width + col0] : 0;
		float nElement1 = ((ty < length) && (col1 < width))? matrix2[ty * width + col1] : 0;
		float nElement2 = ((ty < length) && (col2 < width))? matrix2[ty * width + col2] : 0;
		float nElement3 = ((ty < length) && (col3 < width))? matrix2[ty * width + col3] : 0;

		for (int m = 1; m <= (length + 7) / 8; ++m) {
			mArray[ty][tx] = mElement;
			nArray0[ty][tx] = nElement0;
			nArray1[ty][tx] = nElement1;
			nArray2[ty][tx] = nElement2;
			nArray3[ty][tx] = nElement3;
			__syncthreads();

			mElement = (((m * 8 + tx) < length) && (row < height))? matrix1[row * length + m * 8 + tx] : 0;
			nElement0 = (((m * 8 + ty) < length) && (col0 < width))? matrix2[(m * 8 + ty) * width + col0] : 0;
			nElement1 = (((m * 8 + ty) < length) && (col1 < width))? matrix2[(m * 8 + ty) * width + col1] : 0;
			nElement2 = (((m * 8 + ty) < length) && (col2 < width))? matrix2[(m * 8 + ty) * width + col2] : 0;
			nElement3 = (((m * 8 + ty) < length) && (col3 < width))? matrix2[(m * 8 + ty) * width + col3] : 0;
			
			for (int k = 0; k < 8; ++k) {
				value0 += mArray[ty][k] * nArray0[k][tx];
				value1 += mArray[ty][k] * nArray1[k][tx];
				value2 += mArray[ty][k] * nArray2[k][tx];
				value3 += mArray[ty][k] * nArray3[k][tx];
			}
			__syncthreads();
		}
		if (row < height) {
			if (col0 < width)
				result[row * width + col0] = value0;
			if (col1 < width)
				result[row * width + col1] = value1;
			if (col2 < width)
				result[row * width + col2] = value2;
			if (col3 < width)
				result[row * width + col3] = value3;
		}
}

__global__  void multiplyKernel_8_prefetch_granularity_1x4_2LoopUnrolling(
	float *result, const float *matrix1, const float *matrix2, const int height, const int width, const int length) { 
		__shared__ float mArray[8][8]; 
		__shared__ float nArray0[8][8];
		__shared__ float nArray1[8][8];
		__shared__ float nArray2[8][8];
		__shared__ float nArray3[8][8];
		
		int bx = blockIdx.x; 
		int by = blockIdx.y; 
		int tx = threadIdx.x; 
		int ty = threadIdx.y; 

		int row = by * 8 + ty;
		int col0 = 4 * bx * 8 + tx;
		int col1 = (4 * bx + 1) * 8 + tx;
		int col2 = (4 * bx + 2) * 8 + tx;
		int col3 = (4 * bx + 3) * 8 + tx;

		float value0 = 0;
		float value1 = 0;
		float value2 = 0;
		float value3 = 0;

		float mElement = ((tx < length) && (row < height))? matrix1[row * length + tx] : 0;
		float nElement0 = ((ty < length) && (col0 < width))? matrix2[ty * width + col0] : 0;
		float nElement1 = ((ty < length) && (col1 < width))? matrix2[ty * width + col1] : 0;
		float nElement2 = ((ty < length) && (col2 < width))? matrix2[ty * width + col2] : 0;
		float nElement3 = ((ty < length) && (col3 < width))? matrix2[ty * width + col3] : 0;

		for (int m = 1; m <= (length + 7) / 8; ++m) {
			mArray[ty][tx] = mElement;
			nArray0[ty][tx] = nElement0;
			nArray1[ty][tx] = nElement1;
			nArray2[ty][tx] = nElement2;
			nArray3[ty][tx] = nElement3;
			__syncthreads();

			mElement = (((m * 8 + tx) < length) && (row < height))? matrix1[row * length + m * 8 + tx] : 0;
			nElement0 = (((m * 8 + ty) < length) && (col0 < width))? matrix2[(m * 8 + ty) * width + col0] : 0;
			nElement1 = (((m * 8 + ty) < length) && (col1 < width))? matrix2[(m * 8 + ty) * width + col1] : 0;
			nElement2 = (((m * 8 + ty) < length) && (col2 < width))? matrix2[(m * 8 + ty) * width + col2] : 0;
			nElement3 = (((m * 8 + ty) < length) && (col3 < width))? matrix2[(m * 8 + ty) * width + col3] : 0;
			
			for (int k = 0; k < 4; ++k) {
				value0 += mArray[ty][2 * k] * nArray0[2 * k][tx] + 
					mArray[ty][2 * k + 1] * nArray0[2 * k + 1][tx]; 
				value1 += mArray[ty][2 * k] * nArray1[2 * k][tx] +
					mArray[ty][2 * k + 1] * nArray1[2 * k + 1][tx]; 
				value2 += mArray[ty][2 * k] * nArray2[2 * k][tx] +
					mArray[ty][2 * k + 1] * nArray2[2 * k + 1][tx]; 
				value3 += mArray[ty][2 * k] * nArray3[2 * k][tx] +
					mArray[ty][2 * k + 1] * nArray3[2 * k + 1][tx]; 
			}
			__syncthreads();
		}
		if (row < height) {
			if (col0 < width)
				result[row * width + col0] = value0;
			if (col1 < width)
				result[row * width + col1] = value1;
			if (col2 < width)
				result[row * width + col2] = value2;
			if (col3 < width)
				result[row * width + col3] = value3;
		}
}

__global__  void multiplyKernel_8_prefetch_granularity_1x4_4LoopUnrolling(
	float *result, const float *matrix1, const float *matrix2, const int height, const int width, const int length) { 
		__shared__ float mArray[8][8]; 
		__shared__ float nArray0[8][8];
		__shared__ float nArray1[8][8];
		__shared__ float nArray2[8][8];
		__shared__ float nArray3[8][8];
		
		int bx = blockIdx.x; 
		int by = blockIdx.y; 
		int tx = threadIdx.x; 
		int ty = threadIdx.y; 

		int row = by * 8 + ty;
		int col0 = 4 * bx * 8 + tx;
		int col1 = (4 * bx + 1) * 8 + tx;
		int col2 = (4 * bx + 2) * 8 + tx;
		int col3 = (4 * bx + 3) * 8 + tx;

		float value0 = 0;
		float value1 = 0;
		float value2 = 0;
		float value3 = 0;

		float mElement = ((tx < length) && (row < height))? matrix1[row * length + tx] : 0;
		float nElement0 = ((ty < length) && (col0 < width))? matrix2[ty * width + col0] : 0;
		float nElement1 = ((ty < length) && (col1 < width))? matrix2[ty * width + col1] : 0;
		float nElement2 = ((ty < length) && (col2 < width))? matrix2[ty * width + col2] : 0;
		float nElement3 = ((ty < length) && (col3 < width))? matrix2[ty * width + col3] : 0;

		for (int m = 1; m <= (length + 7) / 8; ++m) {
			mArray[ty][tx] = mElement;
			nArray0[ty][tx] = nElement0;
			nArray1[ty][tx] = nElement1;
			nArray2[ty][tx] = nElement2;
			nArray3[ty][tx] = nElement3;
			__syncthreads();

			mElement = (((m * 8 + tx) < length) && (row < height))? matrix1[row * length + m * 8 + tx] : 0;
			nElement0 = (((m * 8 + ty) < length) && (col0 < width))? matrix2[(m * 8 + ty) * width + col0] : 0;
			nElement1 = (((m * 8 + ty) < length) && (col1 < width))? matrix2[(m * 8 + ty) * width + col1] : 0;
			nElement2 = (((m * 8 + ty) < length) && (col2 < width))? matrix2[(m * 8 + ty) * width + col2] : 0;
			nElement3 = (((m * 8 + ty) < length) && (col3 < width))? matrix2[(m * 8 + ty) * width + col3] : 0;
			
			for (int k = 0; k < 2; ++k) {
				value0 += mArray[ty][4 * k] * nArray0[4 * k][tx] +
					mArray[ty][4 * k + 1] * nArray0[4 * k + 1][tx] + 
					mArray[ty][4 * k + 2] * nArray0[4 * k + 2][tx] + 
					mArray[ty][4 * k + 3] * nArray0[4 * k + 3][tx]; 
				value1 += mArray[ty][4 * k] * nArray1[4 * k][tx] +
					mArray[ty][4 * k + 1] * nArray1[4 * k + 1][tx] + 
					mArray[ty][4 * k + 2] * nArray1[4 * k + 2][tx] + 
					mArray[ty][4 * k + 3] * nArray1[4 * k + 3][tx]; 
				value2 += mArray[ty][4 * k] * nArray2[4 * k][tx] +
					mArray[ty][4 * k + 1] * nArray2[4 * k + 1][tx] + 
					mArray[ty][4 * k + 2] * nArray2[4 * k + 2][tx] + 
					mArray[ty][4 * k + 3] * nArray2[4 * k + 3][tx]; 
				value3 += mArray[ty][4 * k] * nArray3[4 * k][tx] +
					mArray[ty][4 * k + 1] * nArray3[4 * k + 1][tx] + 
					mArray[ty][4 * k + 2] * nArray3[4 * k + 2][tx] + 
					mArray[ty][4 * k + 3] * nArray3[4 * k + 3][tx]; 
			}
			__syncthreads();
		}
		if (row < height) {
			if (col0 < width)
				result[row * width + col0] = value0;
			if (col1 < width)
				result[row * width + col1] = value1;
			if (col2 < width)
				result[row * width + col2] = value2;
			if (col3 < width)
				result[row * width + col3] = value3;
		}
}

__global__  void multiplyKernel_8_prefetch_granularity_1x4_fullLoopUnrolling(
	float *result, const float *matrix1, const float *matrix2, const int height, const int width, const int length) { 
		__shared__ float mArray[8][8]; 
		__shared__ float nArray0[8][8];
		__shared__ float nArray1[8][8];
		__shared__ float nArray2[8][8];
		__shared__ float nArray3[8][8];
		
		int bx = blockIdx.x; 
		int by = blockIdx.y; 
		int tx = threadIdx.x; 
		int ty = threadIdx.y; 

		int row = by * 8 + ty;
		int col0 = 4 * bx * 8 + tx;
		int col1 = (4 * bx + 1) * 8 + tx;
		int col2 = (4 * bx + 2) * 8 + tx;
		int col3 = (4 * bx + 3) * 8 + tx;

		float value0 = 0;
		float value1 = 0;
		float value2 = 0;
		float value3 = 0;

		float mElement = ((tx < length) && (row < height))? matrix1[row * length + tx] : 0;
		float nElement0 = ((ty < length) && (col0 < width))? matrix2[ty * width + col0] : 0;
		float nElement1 = ((ty < length) && (col1 < width))? matrix2[ty * width + col1] : 0;
		float nElement2 = ((ty < length) && (col2 < width))? matrix2[ty * width + col2] : 0;
		float nElement3 = ((ty < length) && (col3 < width))? matrix2[ty * width + col3] : 0;

		for (int m = 1; m <= (length + 7) / 8; ++m) {
			mArray[ty][tx] = mElement;
			nArray0[ty][tx] = nElement0;
			nArray1[ty][tx] = nElement1;
			nArray2[ty][tx] = nElement2;
			nArray3[ty][tx] = nElement3;
			__syncthreads();

			mElement = (((m * 8 + tx) < length) && (row < height))? matrix1[row * length + m * 8 + tx] : 0;
			nElement0 = (((m * 8 + ty) < length) && (col0 < width))? matrix2[(m * 8 + ty) * width + col0] : 0;
			nElement1 = (((m * 8 + ty) < length) && (col1 < width))? matrix2[(m * 8 + ty) * width + col1] : 0;
			nElement2 = (((m * 8 + ty) < length) && (col2 < width))? matrix2[(m * 8 + ty) * width + col2] : 0;
			nElement3 = (((m * 8 + ty) < length) && (col3 < width))? matrix2[(m * 8 + ty) * width + col3] : 0;
			
			value0 += mArray[ty][0] * nArray0[0][tx] +
				mArray[ty][1] * nArray0[1][tx] +
				mArray[ty][2] * nArray0[2][tx] +
				mArray[ty][3] * nArray0[3][tx] +
				mArray[ty][4] * nArray0[4][tx] +
				mArray[ty][5] * nArray0[5][tx] +
				mArray[ty][6] * nArray0[6][tx] +
				mArray[ty][7] * nArray0[7][tx];
			value1 += mArray[ty][0] * nArray1[0][tx] +
				mArray[ty][1] * nArray1[1][tx] +
				mArray[ty][2] * nArray1[2][tx] +
				mArray[ty][3] * nArray1[3][tx] +
				mArray[ty][4] * nArray1[4][tx] +
				mArray[ty][5] * nArray1[5][tx] +
				mArray[ty][6] * nArray1[6][tx] +
				mArray[ty][7] * nArray1[7][tx];
			value2 += mArray[ty][0] * nArray2[0][tx] +
				mArray[ty][1] * nArray2[1][tx] +
				mArray[ty][2] * nArray2[2][tx] +
				mArray[ty][3] * nArray2[3][tx] +
				mArray[ty][4] * nArray2[4][tx] +
				mArray[ty][5] * nArray2[5][tx] +
				mArray[ty][6] * nArray2[6][tx] +
				mArray[ty][7] * nArray2[7][tx];
			value3 += mArray[ty][0] * nArray3[0][tx] +
				mArray[ty][1] * nArray3[1][tx] +
				mArray[ty][2] * nArray3[2][tx] +
				mArray[ty][3] * nArray3[3][tx] +
				mArray[ty][4] * nArray3[4][tx] +
				mArray[ty][5] * nArray3[5][tx] +
				mArray[ty][6] * nArray3[6][tx] +
				mArray[ty][7] * nArray3[7][tx];

			__syncthreads();
		}
		if (row < height) {
			if (col0 < width)
				result[row * width + col0] = value0;
			if (col1 < width)
				result[row * width + col1] = value1;
			if (col2 < width)
				result[row * width + col2] = value2;
			if (col3 < width)
				result[row * width + col3] = value3;
		}
}

__global__  void multiplyKernel_16_noPrefetch_granularity_1x1_noLoopUnrolling(
	float *result, const float *matrix1, const float *matrix2, const int height, const int width, const int length) { 
		__shared__ float mArray[16][16]; 
		__shared__ float nArray[16][16];

		int bx = blockIdx.x; 
		int by = blockIdx.y; 
		int tx = threadIdx.x; 
		int ty = threadIdx.y; 

		int row = by * 16 + ty;
		int col = bx * 16 + tx;

		float value = 0;
		for (int m = 0; m < (length + 15) / 16; ++m) {
			mArray[ty][tx] = (((m * 16 + tx) < length) && (row < height))? matrix1[row * length + m * 16 + tx] : 0;
			nArray[ty][tx] = (((m * 16 + ty) < length) && (col < width))? matrix2[(m * 16 + ty) * width + col] : 0;
			__syncthreads();

			for (int k = 0; k < 16; ++k) {
				value += mArray[ty][k] * nArray[k][tx]; 
			}
			__syncthreads();
		}
		if (row < height && col < width)
			result[row * width + col] = value;
}

__global__  void multiplyKernel_16_noPrefetch_granularity_1x1_2LoopUnrolling(
	float *result, const float *matrix1, const float *matrix2, const int height, const int width, const int length) { 
		__shared__ float mArray[16][16]; 
		__shared__ float nArray[16][16];

		int bx = blockIdx.x; 
		int by = blockIdx.y; 
		int tx = threadIdx.x; 
		int ty = threadIdx.y; 

		int row = by * 16 + ty;
		int col = bx * 16 + tx;

		float value = 0;
		for (int m = 0; m < (length + 15) / 16; ++m) {
			mArray[ty][tx] = (((m * 16 + tx) < length) && (row < height))? matrix1[row * length + m * 16 + tx] : 0;
			nArray[ty][tx] = (((m * 16 + ty) < length) && (col < width))? matrix2[(m * 16 + ty) * width + col] : 0;
			__syncthreads();

			for (int k = 0; k < 8; ++k) {
				value += mArray[ty][2 * k] * nArray[2 * k][tx] + 
					mArray[ty][2 * k + 1] * nArray[2 * k + 1][tx]; 
			}
			__syncthreads();
		}
		if (row < height && col < width)
			result[row * width + col] = value;
}

__global__  void multiplyKernel_16_noPrefetch_granularity_1x1_4LoopUnrolling(
	float *result, const float *matrix1, const float *matrix2, const int height, const int width, const int length) { 
		__shared__ float mArray[16][16]; 
		__shared__ float nArray[16][16];

		int bx = blockIdx.x; 
		int by = blockIdx.y; 
		int tx = threadIdx.x; 
		int ty = threadIdx.y; 

		int row = by * 16 + ty;
		int col = bx * 16 + tx;

		float value = 0;
		for (int m = 0; m < (length + 15) / 16; ++m) {
			mArray[ty][tx] = (((m * 16 + tx) < length) && (row < height))? matrix1[row * length + m * 16 + tx] : 0;
			nArray[ty][tx] = (((m * 16 + ty) < length) && (col < width))? matrix2[(m * 16 + ty) * width + col] : 0;
			__syncthreads();

			for (int k = 0; k < 4; ++k) {
				value += mArray[ty][4 * k] * nArray[4 * k][tx] +
					mArray[ty][4 * k + 1] * nArray[4 * k + 1][tx] +
					mArray[ty][4 * k + 2] * nArray[4 * k + 2][tx] +
					mArray[ty][4 * k + 3] * nArray[4 * k + 3][tx];
			}
			__syncthreads();
		}
		if (row < height && col < width)
			result[row * width + col] = value;
}

__global__  void multiplyKernel_16_noPrefetch_granularity_1x1_fullLoopUnrolling(
	float *result, const float *matrix1, const float *matrix2, const int height, const int width, const int length) { 
		__shared__ float mArray[16][16]; 
		__shared__ float nArray[16][16];

		int bx = blockIdx.x; 
		int by = blockIdx.y; 
		int tx = threadIdx.x; 
		int ty = threadIdx.y; 

		int row = by * 16 + ty;
		int col = bx * 16 + tx;

		float value = 0;
		for (int m = 0; m < (length + 15) / 16; ++m) {
			mArray[ty][tx] = (((m * 16 + tx) < length) && (row < height))? matrix1[row * length + m * 16 + tx] : 0;
			nArray[ty][tx] = (((m * 16 + ty) < length) && (col < width))? matrix2[(m * 16 + ty) * width + col] : 0;
			__syncthreads();

			value += mArray[ty][0] * nArray[0][tx] +
				mArray[ty][1] * nArray[1][tx] +
				mArray[ty][2] * nArray[2][tx] +
				mArray[ty][3] * nArray[3][tx] +
				mArray[ty][4] * nArray[4][tx] +
				mArray[ty][5] * nArray[5][tx] +
				mArray[ty][6] * nArray[6][tx] +
				mArray[ty][7] * nArray[7][tx] +
				mArray[ty][8] * nArray[8][tx] +
				mArray[ty][9] * nArray[9][tx] +
				mArray[ty][10] * nArray[10][tx] +
				mArray[ty][11] * nArray[11][tx] +
				mArray[ty][12] * nArray[12][tx] +
				mArray[ty][13] * nArray[13][tx] +
				mArray[ty][14] * nArray[14][tx] +
				mArray[ty][15] * nArray[15][tx];
			__syncthreads();
		}
		if (row < height && col < width)
			result[row * width + col] = value;
}

__global__  void multiplyKernel_16_noPrefetch_granularity_1x2_noLoopUnrolling(
	float *result, const float *matrix1, const float *matrix2, const int height, const int width, const int length) { 
		__shared__ float mArray[16][16]; 
		__shared__ float nArray0[16][16];
		__shared__ float nArray1[16][16];

		int bx = blockIdx.x; 
		int by = blockIdx.y; 
		int tx = threadIdx.x; 
		int ty = threadIdx.y; 

		int row = by * 16 + ty;
		int col0 = 2 * bx * 16 + tx;
		int col1 = (2 * bx + 1) * 16 + tx;

		float value0 = 0;
		float value1 = 0;
		for (int m = 0; m < (length + 15) / 16; ++m) {
			mArray[ty][tx] = (((m * 16 + tx) < length) && (row < height))? matrix1[row * length + m * 16 + tx] : 0;
			nArray0[ty][tx] = (((m * 16 + ty) < length) && (col0 < width))? matrix2[(m * 16 + ty) * width + col0] : 0;
			nArray1[ty][tx] = (((m * 16 + ty) < length) && (col1 < width))? matrix2[(m * 16 + ty) * width + col1] : 0;
			__syncthreads();

			for (int k = 0; k < 16; ++k) {
				value0 += mArray[ty][k] * nArray0[k][tx];
				value1 += mArray[ty][k] * nArray1[k][tx];
			}
			__syncthreads();
		}
		if (row < height) {
			if (col0 < width)
				result[row * width + col0] = value0;
			if (col1 < width)
				result[row * width + col1] = value1;
		}
}

__global__  void multiplyKernel_16_noPrefetch_granularity_1x2_2LoopUnrolling(
	float *result, const float *matrix1, const float *matrix2, const int height, const int width, const int length) { 
			__shared__ float mArray[16][16]; 
			__shared__ float nArray0[16][16];
			__shared__ float nArray1[16][16];

			int bx = blockIdx.x; 
			int by = blockIdx.y; 
			int tx = threadIdx.x; 
			int ty = threadIdx.y; 

			int row = by * 16 + ty;
			int col0 = 2 * bx * 16 + tx;
			int col1 = (2 * bx + 1) * 16 + tx;

			float value0 = 0;
			float value1 = 0;
			for (int m = 0; m < (length + 15) / 16; ++m) {
				mArray[ty][tx] = (((m * 16 + tx) < length) && (row < height))? matrix1[row * length + m * 16 + tx] : 0;
				nArray0[ty][tx] = (((m * 16 + ty) < length) && (col0 < width))? matrix2[(m * 16 + ty) * width + col0] : 0;
				nArray1[ty][tx] = (((m * 16 + ty) < length) && (col1 < width))? matrix2[(m * 16 + ty) * width + col1] : 0;
				__syncthreads();

				for (int k = 0; k < 8; ++k) {
					value0 += mArray[ty][2 * k] * nArray0[2 * k][tx] + mArray[ty][2 * k + 1] * nArray0[2 * k + 1][tx];
					value1 += mArray[ty][2 * k] * nArray1[2 * k][tx] + mArray[ty][2 * k + 1] * nArray1[2 * k + 1][tx]; 
				}
				__syncthreads();
			}
			if (row < height) {
				if (col0 < width)
					result[row * width + col0] = value0;
				if (col1 < width)
					result[row * width + col1] = value1;
			}
}

__global__  void multiplyKernel_16_noPrefetch_granularity_1x2_4LoopUnrolling(
	float *result, const float *matrix1, const float *matrix2, const int height, const int width, const int length) { 
		__shared__ float mArray[16][16]; 
		__shared__ float nArray0[16][16];
		__shared__ float nArray1[16][16];

		int bx = blockIdx.x; 
		int by = blockIdx.y; 
		int tx = threadIdx.x; 
		int ty = threadIdx.y; 

		int row = by * 16 + ty;
		int col0 = 2 * bx * 16 + tx;
		int col1 = (2 * bx + 1) * 16 + tx;

		float value0 = 0;
		float value1 = 0;
		for (int m = 0; m < (length + 15) / 16; ++m) {
			mArray[ty][tx] = (((m * 16 + tx) < length) && (row < height))? matrix1[row * length + m * 16 + tx] : 0;
			nArray0[ty][tx] = (((m * 16 + ty) < length) && (col0 < width))? matrix2[(m * 16 + ty) * width + col0] : 0;
			nArray1[ty][tx] = (((m * 16 + ty) < length) && (col1 < width))? matrix2[(m * 16 + ty) * width + col1] : 0;
			__syncthreads();

			for (int k = 0; k < 4; ++k) {
				value0 += mArray[ty][4 * k] * nArray0[4 * k][tx] +
					mArray[ty][4 * k + 1] * nArray0[4 * k + 1][tx] +
					mArray[ty][4 * k + 2] * nArray0[4 * k + 2][tx] +
					mArray[ty][4 * k + 3] * nArray0[4 * k + 3][tx];
				value1 += mArray[ty][4 * k] * nArray1[4 * k][tx] +
					mArray[ty][4 * k + 1] * nArray1[4 * k + 1][tx] +
					mArray[ty][4 * k + 2] * nArray1[4 * k + 2][tx] +
					mArray[ty][4 * k + 3] * nArray1[4 * k + 3][tx];
			}
			__syncthreads();
		}
		if (row < height) {
			if (col0 < width)
				result[row * width + col0] = value0;
			if (col1 < width)
				result[row * width + col1] = value1;
		}
}

__global__  void multiplyKernel_16_noPrefetch_granularity_1x2_fullLoopUnrolling(
	float *result, const float *matrix1, const float *matrix2, const int height, const int width, const int length) { 
		__shared__ float mArray[16][16]; 
		__shared__ float nArray0[16][16];
		__shared__ float nArray1[16][16];

		int bx = blockIdx.x; 
		int by = blockIdx.y; 
		int tx = threadIdx.x; 
		int ty = threadIdx.y; 

		int row = by * 16 + ty;
		int col0 = 2 * bx * 16 + tx;
		int col1 = (2 * bx + 1) * 16 + tx;

		float value0 = 0;
		float value1 = 0;
		for (int m = 0; m < (length + 15) / 16; ++m) {
			mArray[ty][tx] = (((m * 16 + tx) < length) && (row < height))? matrix1[row * length + m * 16 + tx] : 0;
			nArray0[ty][tx] = (((m * 16 + ty) < length) && (col0 < width))? matrix2[(m * 16 + ty) * width + col0] : 0;
			nArray1[ty][tx] = (((m * 16 + ty) < length) && (col1 < width))? matrix2[(m * 16 + ty) * width + col1] : 0;
			__syncthreads();

			value0 += mArray[ty][0] * nArray0[0][tx] +
				mArray[ty][1] * nArray0[1][tx] +
				mArray[ty][2] * nArray0[2][tx] +
				mArray[ty][3] * nArray0[3][tx] +
				mArray[ty][4] * nArray0[4][tx] +
				mArray[ty][5] * nArray0[5][tx] +
				mArray[ty][6] * nArray0[6][tx] +
				mArray[ty][7] * nArray0[7][tx] +
				mArray[ty][8] * nArray0[8][tx] +
				mArray[ty][9] * nArray0[9][tx] +
				mArray[ty][10] * nArray0[10][tx] +
				mArray[ty][11] * nArray0[11][tx] +
				mArray[ty][12] * nArray0[12][tx] +
				mArray[ty][13] * nArray0[13][tx] +
				mArray[ty][14] * nArray0[14][tx] +
				mArray[ty][15] * nArray0[15][tx];	
			value1 += mArray[ty][0] * nArray1[0][tx] +
				mArray[ty][1] * nArray1[1][tx] +
				mArray[ty][2] * nArray1[2][tx] +
				mArray[ty][3] * nArray1[3][tx] +
				mArray[ty][4] * nArray1[4][tx] +
				mArray[ty][5] * nArray1[5][tx] +
				mArray[ty][6] * nArray1[6][tx] +
				mArray[ty][7] * nArray1[7][tx] +
				mArray[ty][8] * nArray1[8][tx] +
				mArray[ty][9] * nArray1[9][tx] +
				mArray[ty][10] * nArray1[10][tx] +
				mArray[ty][11] * nArray1[11][tx] +
				mArray[ty][12] * nArray1[12][tx] +
				mArray[ty][13] * nArray1[13][tx] +
				mArray[ty][14] * nArray1[14][tx] +
				mArray[ty][15] * nArray1[15][tx];
			__syncthreads();
		}
		if (row < height) {
			if (col0 < width)
				result[row * width + col0] = value0;
			if (col1 < width)
				result[row * width + col1] = value1;
		}
}

__global__  void multiplyKernel_16_noPrefetch_granularity_1x4_noLoopUnrolling(
	float *result, const float *matrix1, const float *matrix2, const int height, const int width, const int length) { 
		__shared__ float mArray[16][16]; 
		__shared__ float nArray0[16][16];
		__shared__ float nArray1[16][16];
		__shared__ float nArray2[16][16];
		__shared__ float nArray3[16][16];

		int bx = blockIdx.x; 
		int by = blockIdx.y; 
		int tx = threadIdx.x; 
		int ty = threadIdx.y; 

		int row = by * 16 + ty;
		int col0 = 4 * bx * 16 + tx;
		int col1 = (4 * bx + 1) * 16 + tx;
		int col2 = (4 * bx + 2) * 16 + tx;
		int col3 = (4 * bx + 3) * 16 + tx;

		float value0 = 0;
		float value1 = 0;
		float value2 = 0;
		float value3 = 0;
		for (int m = 0; m < (length + 15) / 16; ++m) {
			mArray[ty][tx] = (((m * 16 + tx) < length) && (row < height))? matrix1[row * length + m * 16 + tx] : 0;
			nArray0[ty][tx] = (((m * 16 + ty) < length) && (col0 < width))? matrix2[(m * 16 + ty) * width + col0] : 0;
			nArray1[ty][tx] = (((m * 16 + ty) < length) && (col1 < width))? matrix2[(m * 16 + ty) * width + col1] : 0;
			nArray2[ty][tx] = (((m * 16 + ty) < length) && (col2 < width))? matrix2[(m * 16 + ty) * width + col2] : 0;
			nArray3[ty][tx] = (((m * 16 + ty) < length) && (col3 < width))? matrix2[(m * 16 + ty) * width + col3] : 0;
			__syncthreads();

			for (int k = 0; k < 16; ++k) {
				value0 += mArray[ty][k] * nArray0[k][tx];
				value1 += mArray[ty][k] * nArray1[k][tx];
				value2 += mArray[ty][k] * nArray2[k][tx];
				value3 += mArray[ty][k] * nArray3[k][tx];
			}
			__syncthreads();
		}
		if (row < height) {
			if (col0 < width)
				result[row * width + col0] = value0;
			if (col1 < width)
				result[row * width + col1] = value1;
			if (col2 < width)
				result[row * width + col2] = value2;
			if (col3 < width)
				result[row * width + col3] = value3;
		}
}

__global__  void multiplyKernel_16_noPrefetch_granularity_1x4_2LoopUnrolling(
	float *result, const float *matrix1, const float *matrix2, const int height, const int width, const int length) { 
		__shared__ float mArray[16][16]; 
		__shared__ float nArray0[16][16];
		__shared__ float nArray1[16][16];
		__shared__ float nArray2[16][16];
		__shared__ float nArray3[16][16];

		int bx = blockIdx.x; 
		int by = blockIdx.y; 
		int tx = threadIdx.x; 
		int ty = threadIdx.y; 

		int row = by * 16 + ty;
		int col0 = 4 * bx * 16 + tx;
		int col1 = (4 * bx + 1) * 16 + tx;
		int col2 = (4 * bx + 2) * 16 + tx;
		int col3 = (4 * bx + 3) * 16 + tx;

		float value0 = 0;
		float value1 = 0;
		float value2 = 0;
		float value3 = 0;
		for (int m = 0; m < (length + 15) / 16; ++m) {
			mArray[ty][tx] = (((m * 16 + tx) < length) && (row < height))? matrix1[row * length + m * 16 + tx] : 0;
			nArray0[ty][tx] = (((m * 16 + ty) < length) && (col0 < width))? matrix2[(m * 16 + ty) * width + col0] : 0;
			nArray1[ty][tx] = (((m * 16 + ty) < length) && (col1 < width))? matrix2[(m * 16 + ty) * width + col1] : 0;
			nArray2[ty][tx] = (((m * 16 + ty) < length) && (col2 < width))? matrix2[(m * 16 + ty) * width + col2] : 0;
			nArray3[ty][tx] = (((m * 16 + ty) < length) && (col3 < width))? matrix2[(m * 16 + ty) * width + col3] : 0;
			__syncthreads();

			for (int k = 0; k < 8; ++k) {
				value0 += mArray[ty][2 * k] * nArray0[2 * k][tx] +
					mArray[ty][2 * k + 1] * nArray0[2 * k + 1][tx];
				value1 += mArray[ty][2 * k] * nArray1[2 * k][tx] +
					mArray[ty][2 * k + 1] * nArray1[2 * k + 1][tx];
				value2 += mArray[ty][2 * k] * nArray2[2 * k][tx] +
					mArray[ty][2 * k + 1] * nArray2[2 * k + 1][tx];
				value3 += mArray[ty][2 * k] * nArray3[2 * k][tx] +
					mArray[ty][2 * k + 1] * nArray3[2 * k + 1][tx];
			}
			__syncthreads();
		}
		if (row < height) {
			if (col0 < width)
				result[row * width + col0] = value0;
			if (col1 < width)
				result[row * width + col1] = value1;
			if (col2 < width)
				result[row * width + col2] = value2;
			if (col3 < width)
				result[row * width + col3] = value3;
		}
}

__global__  void multiplyKernel_16_noPrefetch_granularity_1x4_4LoopUnrolling(
	float *result, const float *matrix1, const float *matrix2, const int height, const int width, const int length) { 
		__shared__ float mArray[16][16]; 
		__shared__ float nArray0[16][16];
		__shared__ float nArray1[16][16];
		__shared__ float nArray2[16][16];
		__shared__ float nArray3[16][16];

		int bx = blockIdx.x; 
		int by = blockIdx.y; 
		int tx = threadIdx.x; 
		int ty = threadIdx.y; 

		int row = by * 16 + ty;
		int col0 = 4 * bx * 16 + tx;
		int col1 = (4 * bx + 1) * 16 + tx;
		int col2 = (4 * bx + 2) * 16 + tx;
		int col3 = (4 * bx + 3) * 16 + tx;

		float value0 = 0;
		float value1 = 0;
		float value2 = 0;
		float value3 = 0;
		for (int m = 0; m < (length + 15) / 16; ++m) {
			mArray[ty][tx] = (((m * 16 + tx) < length) && (row < height))? matrix1[row * length + m * 16 + tx] : 0;
			nArray0[ty][tx] = (((m * 16 + ty) < length) && (col0 < width))? matrix2[(m * 16 + ty) * width + col0] : 0;
			nArray1[ty][tx] = (((m * 16 + ty) < length) && (col1 < width))? matrix2[(m * 16 + ty) * width + col1] : 0;
			nArray2[ty][tx] = (((m * 16 + ty) < length) && (col2 < width))? matrix2[(m * 16 + ty) * width + col2] : 0;
			nArray3[ty][tx] = (((m * 16 + ty) < length) && (col3 < width))? matrix2[(m * 16 + ty) * width + col3] : 0;
			__syncthreads();

			for (int k = 0; k < 4; ++k) {
				value0 += mArray[ty][4 * k] * nArray0[4 * k][tx] +
					mArray[ty][4 * k + 1] * nArray0[4 * k + 1][tx] +
					mArray[ty][4 * k + 2] * nArray0[4 * k + 2][tx] +
					mArray[ty][4 * k + 3] * nArray0[4 * k + 3][tx];
				value1 += mArray[ty][4 * k] * nArray1[4 * k][tx] +
					mArray[ty][4 * k + 1] * nArray1[4 * k + 1][tx] +
					mArray[ty][4 * k + 2] * nArray1[4 * k + 2][tx] +
					mArray[ty][4 * k + 3] * nArray1[4 * k + 3][tx];
				value2 += mArray[ty][4 * k] * nArray2[4 * k][tx] +
					mArray[ty][4 * k + 1] * nArray2[4 * k + 1][tx] +
					mArray[ty][4 * k + 2] * nArray2[4 * k + 2][tx] +
					mArray[ty][4 * k + 3] * nArray2[4 * k + 3][tx];
				value3 += mArray[ty][4 * k] * nArray3[4 * k][tx] +
					mArray[ty][4 * k + 1] * nArray3[4 * k + 1][tx] +
					mArray[ty][4 * k + 2] * nArray3[4 * k + 2][tx] +
					mArray[ty][4 * k + 3] * nArray3[4 * k + 3][tx];
			}
			__syncthreads();
		}
		if (row < height) {
			if (col0 < width)
				result[row * width + col0] = value0;
			if (col1 < width)
				result[row * width + col1] = value1;
			if (col2 < width)
				result[row * width + col2] = value2;
			if (col3 < width)
				result[row * width + col3] = value3;
		}
}

__global__  void multiplyKernel_16_noPrefetch_granularity_1x4_fullLoopUnrolling(
	float *result, const float *matrix1, const float *matrix2, const int height, const int width, const int length) { 
		__shared__ float mArray[16][16]; 
		__shared__ float nArray0[16][16];
		__shared__ float nArray1[16][16];
		__shared__ float nArray2[16][16];
		__shared__ float nArray3[16][16];

		int bx = blockIdx.x; 
		int by = blockIdx.y; 
		int tx = threadIdx.x; 
		int ty = threadIdx.y; 

		int row = by * 16 + ty;
		int col0 = 4 * bx * 16 + tx;
		int col1 = (4 * bx + 1) * 16 + tx;
		int col2 = (4 * bx + 2) * 16 + tx;
		int col3 = (4 * bx + 3) * 16 + tx;

		float value0 = 0;
		float value1 = 0;
		float value2 = 0;
		float value3 = 0;
		for (int m = 0; m < (length + 15) / 16; ++m) {
			mArray[ty][tx] = (((m * 16 + tx) < length) && (row < height))? matrix1[row * length + m * 16 + tx] : 0;
			nArray0[ty][tx] = (((m * 16 + ty) < length) && (col0 < width))? matrix2[(m * 16 + ty) * width + col0] : 0;
			nArray1[ty][tx] = (((m * 16 + ty) < length) && (col1 < width))? matrix2[(m * 16 + ty) * width + col1] : 0;
			nArray2[ty][tx] = (((m * 16 + ty) < length) && (col2 < width))? matrix2[(m * 16 + ty) * width + col2] : 0;
			nArray3[ty][tx] = (((m * 16 + ty) < length) && (col3 < width))? matrix2[(m * 16 + ty) * width + col3] : 0;
			__syncthreads();

			value0 += mArray[ty][0] * nArray0[0][tx] +
				mArray[ty][1] * nArray0[1][tx] +
				mArray[ty][2] * nArray0[2][tx] +
				mArray[ty][3] * nArray0[3][tx] +
				mArray[ty][4] * nArray0[4][tx] +
				mArray[ty][5] * nArray0[5][tx] +
				mArray[ty][6] * nArray0[6][tx] +
				mArray[ty][7] * nArray0[7][tx] +
				mArray[ty][8] * nArray0[8][tx] +
				mArray[ty][9] * nArray0[9][tx] +
				mArray[ty][10] * nArray0[10][tx] +
				mArray[ty][11] * nArray0[11][tx] +
				mArray[ty][12] * nArray0[12][tx] +
				mArray[ty][13] * nArray0[13][tx] +
				mArray[ty][14] * nArray0[14][tx] +
				mArray[ty][15] * nArray0[15][tx];	
			value1 += mArray[ty][0] * nArray1[0][tx] +
				mArray[ty][1] * nArray1[1][tx] +
				mArray[ty][2] * nArray1[2][tx] +
				mArray[ty][3] * nArray1[3][tx] +
				mArray[ty][4] * nArray1[4][tx] +
				mArray[ty][5] * nArray1[5][tx] +
				mArray[ty][6] * nArray1[6][tx] +
				mArray[ty][7] * nArray1[7][tx] +
				mArray[ty][8] * nArray1[8][tx] +
				mArray[ty][9] * nArray1[9][tx] +
				mArray[ty][10] * nArray1[10][tx] +
				mArray[ty][11] * nArray1[11][tx] +
				mArray[ty][12] * nArray1[12][tx] +
				mArray[ty][13] * nArray1[13][tx] +
				mArray[ty][14] * nArray1[14][tx] +
				mArray[ty][15] * nArray1[15][tx];
			value2 += mArray[ty][0] * nArray2[0][tx] +
				mArray[ty][1] * nArray2[1][tx] +
				mArray[ty][2] * nArray2[2][tx] +
				mArray[ty][3] * nArray2[3][tx] +
				mArray[ty][4] * nArray2[4][tx] +
				mArray[ty][5] * nArray2[5][tx] +
				mArray[ty][6] * nArray2[6][tx] +
				mArray[ty][7] * nArray2[7][tx] +
				mArray[ty][8] * nArray2[8][tx] +
				mArray[ty][9] * nArray2[9][tx] +
				mArray[ty][10] * nArray2[10][tx] +
				mArray[ty][11] * nArray2[11][tx] +
				mArray[ty][12] * nArray2[12][tx] +
				mArray[ty][13] * nArray2[13][tx] +
				mArray[ty][14] * nArray2[14][tx] +
				mArray[ty][15] * nArray2[15][tx];
			value3 += mArray[ty][0] * nArray3[0][tx] +
				mArray[ty][1] * nArray3[1][tx] +
				mArray[ty][2] * nArray3[2][tx] +
				mArray[ty][3] * nArray3[3][tx] +
				mArray[ty][4] * nArray3[4][tx] +
				mArray[ty][5] * nArray3[5][tx] +
				mArray[ty][6] * nArray3[6][tx] +
				mArray[ty][7] * nArray3[7][tx] +
				mArray[ty][8] * nArray3[8][tx] +
				mArray[ty][9] * nArray3[9][tx] +
				mArray[ty][10] * nArray3[10][tx] +
				mArray[ty][11] * nArray3[11][tx] +
				mArray[ty][12] * nArray3[12][tx] +
				mArray[ty][13] * nArray3[13][tx] +
				mArray[ty][14] * nArray3[14][tx] +
				mArray[ty][15] * nArray3[15][tx];
			__syncthreads();
		}
		if (row < height) {
			if (col0 < width)
				result[row * width + col0] = value0;
			if (col1 < width)
				result[row * width + col1] = value1;
			if (col2 < width)
				result[row * width + col2] = value2;
			if (col3 < width)
				result[row * width + col3] = value3;
		}
}

__global__  void multiplyKernel_16_prefetch_granularity_1x1_noLoopUnrolling(
	float *result, const float *matrix1, const float *matrix2, const int height, const int width, const int length) { 
			__shared__ float mArray[16][16]; 
			__shared__ float nArray[16][16];

			int bx = blockIdx.x; 
			int by = blockIdx.y; 
			int tx = threadIdx.x; 
			int ty = threadIdx.y; 

			int row = by * 16 + ty;
			int col = bx * 16 + tx;

			float value = 0;
			float mElement = ((tx < length) && (row < height))? matrix1[row * length + tx] : 0;
			float nElement = ((ty < length) && (col < width))? matrix2[ty * width + col] : 0;
			for (int m = 1; m <= (length + 15) / 16; ++m) {
				mArray[ty][tx] = mElement;
				nArray[ty][tx] = nElement;
				__syncthreads();

				mElement = (((m * 16 + tx) < length) && (row < height))? matrix1[row * length + m * 16 + tx] : 0;
				nElement = (((m * 16 + ty) < length) && (col < width))? matrix2[(m * 16 + ty) * width + col] : 0;
				for (int k = 0; k < 16; ++k) {
					value += mArray[ty][k] * nArray[k][tx]; 
				}
				__syncthreads();
			}
			if (row < height && col < width)
				result[row * width + col] = value;
}

__global__  void multiplyKernel_16_prefetch_granularity_1x1_2LoopUnrolling(
	float *result, const float *matrix1, const float *matrix2, const int height, const int width, const int length) { 
		__shared__ float mArray[16][16]; 
		__shared__ float nArray[16][16];

		int bx = blockIdx.x; 
		int by = blockIdx.y; 
		int tx = threadIdx.x; 
		int ty = threadIdx.y; 

		int row = by * 16 + ty;
		int col = bx * 16 + tx;

		float value = 0;
		float mElement = ((tx < length) && (row < height))? matrix1[row * length + tx] : 0;
		float nElement = ((ty < length) && (col < width))? matrix2[ty * width + col] : 0;
		for (int m = 1; m <= (length + 15) / 16; ++m) {
			mArray[ty][tx] = mElement;
			nArray[ty][tx] = nElement;
			__syncthreads();

			mElement = (((m * 16 + tx) < length) && (row < height))? matrix1[row * length + m * 16 + tx] : 0;
			nElement = (((m * 16 + ty) < length) && (col < width))? matrix2[(m * 16 + ty) * width + col] : 0;
			for (int k = 0; k < 8; ++k) {
				value += mArray[ty][2 * k] * nArray[2 * k][tx] + 
					mArray[ty][2 * k + 1] * nArray[2 * k + 1][tx]; 
			}
			__syncthreads();
		}
		if (row < height && col < width)
			result[row * width + col] = value;
}

__global__  void multiplyKernel_16_prefetch_granularity_1x1_4LoopUnrolling(
	float *result, const float *matrix1, const float *matrix2, const int height, const int width, const int length) { 
		__shared__ float mArray[16][16]; 
		__shared__ float nArray[16][16];

		int bx = blockIdx.x; 
		int by = blockIdx.y; 
		int tx = threadIdx.x; 
		int ty = threadIdx.y; 

		int row = by * 16 + ty;
		int col = bx * 16 + tx;

		float value = 0;
		float mElement = ((tx < length) && (row < height))? matrix1[row * length + tx] : 0;
		float nElement = ((ty < length) && (col < width))? matrix2[ty * width + col] : 0;
		for (int m = 1; m <= (length + 15) / 16; ++m) {
			mArray[ty][tx] = mElement;
			nArray[ty][tx] = nElement;
			__syncthreads();

			mElement = (((m * 16 + tx) < length) && (row < height))? matrix1[row * length + m * 16 + tx] : 0;
			nElement = (((m * 16 + ty) < length) && (col < width))? matrix2[(m * 16 + ty) * width + col] : 0;
			for (int k = 0; k < 4; ++k) {
				value += mArray[ty][4 * k] * nArray[4 * k][tx] +
					mArray[ty][4 * k + 1] * nArray[4 * k + 1][tx] +
					mArray[ty][4 * k + 2] * nArray[4 * k + 2][tx] +
					mArray[ty][4 * k + 3] * nArray[4 * k + 3][tx];
			}
			__syncthreads();
		}
		if (row < height && col < width)
			result[row * width + col] = value;
}

__global__  void multiplyKernel_16_prefetch_granularity_1x1_fullLoopUnrolling(
	float *result, const float *matrix1, const float *matrix2, const int height, const int width, const int length) { 
		__shared__ float mArray[16][16]; 
		__shared__ float nArray[16][16];

		int bx = blockIdx.x; 
		int by = blockIdx.y; 
		int tx = threadIdx.x; 
		int ty = threadIdx.y; 

		int row = by * 16 + ty;
		int col = bx * 16 + tx;

		float value = 0;
		float mElement = ((tx < length) && (row < height))? matrix1[row * length + tx] : 0;
		float nElement = ((ty < length) && (col < width))? matrix2[ty * width + col] : 0;
		for (int m = 1; m <= (length + 15) / 16; ++m) {
			mArray[ty][tx] = mElement;
			nArray[ty][tx] = nElement;
			__syncthreads();

			mElement = (((m * 16 + tx) < length) && (row < height))? matrix1[row * length + m * 16 + tx] : 0;
			nElement = (((m * 16 + ty) < length) && (col < width))? matrix2[(m * 16 + ty) * width + col] : 0;
			value += mArray[ty][0] * nArray[0][tx] +
				mArray[ty][1] * nArray[1][tx] +
				mArray[ty][2] * nArray[2][tx] +
				mArray[ty][3] * nArray[3][tx] +
				mArray[ty][4] * nArray[4][tx] +
				mArray[ty][5] * nArray[5][tx] +
				mArray[ty][6] * nArray[6][tx] +
				mArray[ty][7] * nArray[7][tx] +
				mArray[ty][8] * nArray[8][tx] +
				mArray[ty][9] * nArray[9][tx] +
				mArray[ty][10] * nArray[10][tx] +
				mArray[ty][11] * nArray[11][tx] +
				mArray[ty][12] * nArray[12][tx] +
				mArray[ty][13] * nArray[13][tx] +
				mArray[ty][14] * nArray[14][tx] +
				mArray[ty][15] * nArray[15][tx];
			__syncthreads();
		}
		if (row < height && col < width)
			result[row * width + col] = value;
}

__global__  void multiplyKernel_16_prefetch_granularity_1x2_noLoopUnrolling(
	float *result, const float *matrix1, const float *matrix2, const int height, const int width, const int length) { 
		__shared__ float mArray[16][16]; 
		__shared__ float nArray0[16][16];
		__shared__ float nArray1[16][16];

		int bx = blockIdx.x; 
		int by = blockIdx.y; 
		int tx = threadIdx.x; 
		int ty = threadIdx.y; 

		int row = by * 16 + ty;
		int col0 = 2 * bx * 16 + tx;
		int col1 = (2 * bx + 1) * 16 + tx;

		float value0 = 0;
		float value1 = 0;

		float mElement = ((tx < length) && (row < height))? matrix1[row * length + tx] : 0;
		float nElement0 = ((ty < length) && (col0 < width))? matrix2[ty * width + col0] : 0;
		float nElement1 = ((ty < length) && (col1 < width))? matrix2[ty * width + col1] : 0;

		for (int m = 1; m <= (length + 15) / 16; ++m) {
			mArray[ty][tx] = mElement;
			nArray0[ty][tx] = nElement0;
			nArray1[ty][tx] = nElement1;
			__syncthreads();

			mElement = (((m * 16 + tx) < length) && (row < height))? matrix1[row * length + m * 16 + tx] : 0;
			nElement0 = (((m * 16 + ty) < length) && (col0 < width))? matrix2[(m * 16 + ty) * width + col0] : 0;
			nElement1 = (((m * 16 + ty) < length) && (col1 < width))? matrix2[(m * 16 + ty) * width + col1] : 0;
		
			for (int k = 0; k < 16; ++k) {
				value0 += mArray[ty][k] * nArray0[k][tx];
				value1 += mArray[ty][k] * nArray1[k][tx];
			}
			__syncthreads();
		}
		if (row < height) {
			if (col0 < width)
				result[row * width + col0] = value0;
			if (col1 < width)
				result[row * width + col1] = value1;
		}
}

__global__  void multiplyKernel_16_prefetch_granularity_1x2_2LoopUnrolling(
	float *result, const float *matrix1, const float *matrix2, const int height, const int width, const int length) { 
		__shared__ float mArray[16][16]; 
		__shared__ float nArray0[16][16];
		__shared__ float nArray1[16][16];

		int bx = blockIdx.x; 
		int by = blockIdx.y; 
		int tx = threadIdx.x; 
		int ty = threadIdx.y; 

		int row = by * 16 + ty;
		int col0 = 2 * bx * 16 + tx;
		int col1 = (2 * bx + 1) * 16 + tx;

		float value0 = 0;
		float value1 = 0;

		float mElement = ((tx < length) && (row < height))? matrix1[row * length + tx] : 0;
		float nElement0 = ((ty < length) && (col0 < width))? matrix2[ty * width + col0] : 0;
		float nElement1 = ((ty < length) && (col1 < width))? matrix2[ty * width + col1] : 0;

		for (int m = 1; m <= (length + 15) / 16; ++m) {
			mArray[ty][tx] = mElement;
			nArray0[ty][tx] = nElement0;
			nArray1[ty][tx] = nElement1;
			__syncthreads();
			
			mElement = (((m * 16 + tx) < length) && (row < height))? matrix1[row * length + m * 16 + tx] : 0;
			nElement0 = (((m * 16 + ty) < length) && (col0 < width))? matrix2[(m * 16 + ty) * width + col0] : 0;
			nElement1 = (((m * 16 + ty) < length) && (col1 < width))? matrix2[(m * 16 + ty) * width + col1] : 0;
		
			for (int k = 0; k < 8; ++k) {
				value0 += mArray[ty][2 * k] * nArray0[2 * k][tx] + 
					mArray[ty][2 * k + 1] * nArray0[2 * k + 1][tx];
				value1 += mArray[ty][2 * k] * nArray1[2 * k][tx] + 
					mArray[ty][2 * k + 1] * nArray1[2 * k + 1][tx]; 
			}
			__syncthreads();
		}
		if (row < height) {
			if (col0 < width)
				result[row * width + col0] = value0;
			if (col1 < width)
				result[row * width + col1] = value1;
		}
}

__global__  void multiplyKernel_16_prefetch_granularity_1x2_4LoopUnrolling(
	float *result, const float *matrix1, const float *matrix2, const int height, const int width, const int length) { 
		__shared__ float mArray[16][16]; 
		__shared__ float nArray0[16][16];
		__shared__ float nArray1[16][16];

		int bx = blockIdx.x; 
		int by = blockIdx.y; 
		int tx = threadIdx.x; 
		int ty = threadIdx.y; 

		int row = by * 16 + ty;
		int col0 = 2 * bx * 16 + tx;
		int col1 = (2 * bx + 1) * 16 + tx;

		float value0 = 0;
		float value1 = 0;

		float mElement = ((tx < length) && (row < height))? matrix1[row * length + tx] : 0;
		float nElement0 = ((ty < length) && (col0 < width))? matrix2[ty * width + col0] : 0;
		float nElement1 = ((ty < length) && (col1 < width))? matrix2[ty * width + col1] : 0;

		for (int m = 1; m <= (length + 15) / 16; ++m) {
			mArray[ty][tx] = mElement;
			nArray0[ty][tx] = nElement0;
			nArray1[ty][tx] = nElement1;
			__syncthreads();

			mElement = (((m * 16 + tx) < length) && (row < height))? matrix1[row * length + m * 16 + tx] : 0;
			nElement0 = (((m * 16 + ty) < length) && (col0 < width))? matrix2[(m * 16 + ty) * width + col0] : 0;
			nElement1 = (((m * 16 + ty) < length) && (col1 < width))? matrix2[(m * 16 + ty) * width + col1] : 0;
		
			for (int k = 0; k < 4; ++k) {
				value0 += mArray[ty][4 * k] * nArray0[4 * k][tx] +
					mArray[ty][4 * k + 1] * nArray0[4 * k + 1][tx] +
					mArray[ty][4 * k + 2] * nArray0[4 * k + 2][tx] +
					mArray[ty][4 * k + 3] * nArray0[4 * k + 3][tx];
				value1 += mArray[ty][4 * k] * nArray1[4 * k][tx] +
					mArray[ty][4 * k + 1] * nArray1[4 * k + 1][tx] +
					mArray[ty][4 * k + 2] * nArray1[4 * k + 2][tx] +
					mArray[ty][4 * k + 3] * nArray1[4 * k + 3][tx];
			}
			__syncthreads();
		}
		if (row < height) {
			if (col0 < width)
				result[row * width + col0] = value0;
			if (col1 < width)
				result[row * width + col1] = value1;
		}
}

__global__  void multiplyKernel_16_prefetch_granularity_1x2_fullLoopUnrolling(
	float *result, const float *matrix1, const float *matrix2, const int height, const int width, const int length) { 
		__shared__ float mArray[16][16]; 
		__shared__ float nArray0[16][16];
		__shared__ float nArray1[16][16];

		int bx = blockIdx.x; 
		int by = blockIdx.y; 
		int tx = threadIdx.x; 
		int ty = threadIdx.y; 

		int row = by * 16 + ty;
		int col0 = 2 * bx * 16 + tx;
		int col1 = (2 * bx + 1) * 16 + tx;

		float value0 = 0;
		float value1 = 0;

		float mElement = ((tx < length) && (row < height))? matrix1[row * length + tx] : 0;
		float nElement0 = ((ty < length) && (col0 < width))? matrix2[ty * width + col0] : 0;
		float nElement1 = ((ty < length) && (col1 < width))? matrix2[ty * width + col1] : 0;
		for (int m = 1; m <= (length + 15) / 16; ++m) {
			mArray[ty][tx] = mElement;
			nArray0[ty][tx] = nElement0;
			nArray1[ty][tx] = nElement1;
			__syncthreads();

			mElement = (((m * 16 + tx) < length) && (row < height))? matrix1[row * length + m * 16 + tx] : 0;
			nElement0 = (((m * 16 + ty) < length) && (col0 < width))? matrix2[(m * 16 + ty) * width + col0] : 0;
			nElement1 = (((m * 16 + ty) < length) && (col1 < width))? matrix2[(m * 16 + ty) * width + col1] : 0;
		
			value0 += mArray[ty][0] * nArray0[0][tx] +
				mArray[ty][1] * nArray0[1][tx] +
				mArray[ty][2] * nArray0[2][tx] +
				mArray[ty][3] * nArray0[3][tx] +
				mArray[ty][4] * nArray0[4][tx] +
				mArray[ty][5] * nArray0[5][tx] +
				mArray[ty][6] * nArray0[6][tx] +
				mArray[ty][7] * nArray0[7][tx] +
				mArray[ty][8] * nArray0[8][tx] +
				mArray[ty][9] * nArray0[9][tx] +
				mArray[ty][10] * nArray0[10][tx] +
				mArray[ty][11] * nArray0[11][tx] +
				mArray[ty][12] * nArray0[12][tx] +
				mArray[ty][13] * nArray0[13][tx] +
				mArray[ty][14] * nArray0[14][tx] +
				mArray[ty][15] * nArray0[15][tx];	
			value1 += mArray[ty][0] * nArray1[0][tx] +
				mArray[ty][1] * nArray1[1][tx] +
				mArray[ty][2] * nArray1[2][tx] +
				mArray[ty][3] * nArray1[3][tx] +
				mArray[ty][4] * nArray1[4][tx] +
				mArray[ty][5] * nArray1[5][tx] +
				mArray[ty][6] * nArray1[6][tx] +
				mArray[ty][7] * nArray1[7][tx] +
				mArray[ty][8] * nArray1[8][tx] +
				mArray[ty][9] * nArray1[9][tx] +
				mArray[ty][10] * nArray1[10][tx] +
				mArray[ty][11] * nArray1[11][tx] +
				mArray[ty][12] * nArray1[12][tx] +
				mArray[ty][13] * nArray1[13][tx] +
				mArray[ty][14] * nArray1[14][tx] +
				mArray[ty][15] * nArray1[15][tx];
			__syncthreads();
		}
		if (row < height) {
			if (col0 < width)
				result[row * width + col0] = value0;
			if (col1 < width)
				result[row * width + col1] = value1;
		}
}

__global__  void multiplyKernel_16_prefetch_granularity_1x4_noLoopUnrolling(
	float *result, const float *matrix1, const float *matrix2, const int height, const int width, const int length) { 
		__shared__ float mArray[16][16]; 
		__shared__ float nArray0[16][16];
		__shared__ float nArray1[16][16];
		__shared__ float nArray2[16][16];
		__shared__ float nArray3[16][16];

		int bx = blockIdx.x; 
		int by = blockIdx.y; 
		int tx = threadIdx.x; 
		int ty = threadIdx.y; 

		int row = by * 16 + ty;
		int col0 = 4 * bx * 16 + tx;
		int col1 = (4 * bx + 1) * 16 + tx;
		int col2 = (4 * bx + 2) * 16 + tx;
		int col3 = (4 * bx + 3) * 16 + tx;

		float value0 = 0;
		float value1 = 0;
		float value2 = 0;
		float value3 = 0;

		float mElement = ((tx < length) && (row < height))? matrix1[row * length + tx] : 0;
		float nElement0 = ((ty < length) && (col0 < width))? matrix2[ty * width + col0] : 0;
		float nElement1 = ((ty < length) && (col1 < width))? matrix2[ty * width + col1] : 0;
		float nElement2 = ((ty < length) && (col2 < width))? matrix2[ty * width + col2] : 0;
		float nElement3 = ((ty < length) && (col3 < width))? matrix2[ty * width + col3] : 0;

		for (int m = 1; m <= (length + 15) / 16; ++m) {
			mArray[ty][tx] = mElement;
			nArray0[ty][tx] = nElement0;
			nArray1[ty][tx] = nElement1;
			nArray2[ty][tx] = nElement2;
			nArray3[ty][tx] = nElement3;
			__syncthreads();

			mElement = (((m * 16 + tx) < length) && (row < height))? matrix1[row * length + m * 16 + tx] : 0;
			nElement0 = (((m * 16 + ty) < length) && (col0 < width))? matrix2[(m * 16 + ty) * width + col0] : 0;
			nElement1 = (((m * 16 + ty) < length) && (col1 < width))? matrix2[(m * 16 + ty) * width + col1] : 0;
			nElement2 = (((m * 16 + ty) < length) && (col2 < width))? matrix2[(m * 16 + ty) * width + col2] : 0;
			nElement3 = (((m * 16 + ty) < length) && (col3 < width))? matrix2[(m * 16 + ty) * width + col3] : 0;
		
			for (int k = 0; k < 16; ++k) {
				value0 += mArray[ty][k] * nArray0[k][tx];
				value1 += mArray[ty][k] * nArray1[k][tx];
				value2 += mArray[ty][k] * nArray2[k][tx];
				value3 += mArray[ty][k] * nArray3[k][tx];
			}
			__syncthreads();
		}
		if (row < height) {
			if (col0 < width)
				result[row * width + col0] = value0;
			if (col1 < width)
				result[row * width + col1] = value1;
			if (col2 < width)
				result[row * width + col2] = value2;
			if (col3 < width)
				result[row * width + col3] = value3;
		}
}

__global__  void multiplyKernel_16_prefetch_granularity_1x4_2LoopUnrolling(
	float *result, const float *matrix1, const float *matrix2, const int height, const int width, const int length) { 
		__shared__ float mArray[16][16]; 
		__shared__ float nArray0[16][16];
		__shared__ float nArray1[16][16];
		__shared__ float nArray2[16][16];
		__shared__ float nArray3[16][16];

		int bx = blockIdx.x; 
		int by = blockIdx.y; 
		int tx = threadIdx.x; 
		int ty = threadIdx.y; 

		int row = by * 16 + ty;
		int col0 = 4 * bx * 16 + tx;
		int col1 = (4 * bx + 1) * 16 + tx;
		int col2 = (4 * bx + 2) * 16 + tx;
		int col3 = (4 * bx + 3) * 16 + tx;

		float value0 = 0;
		float value1 = 0;
		float value2 = 0;
		float value3 = 0;

		float mElement = ((tx < length) && (row < height))? matrix1[row * length + tx] : 0;
		float nElement0 = ((ty < length) && (col0 < width))? matrix2[ty * width + col0] : 0;
		float nElement1 = ((ty < length) && (col1 < width))? matrix2[ty * width + col1] : 0;
		float nElement2 = ((ty < length) && (col2 < width))? matrix2[ty * width + col2] : 0;
		float nElement3 = ((ty < length) && (col3 < width))? matrix2[ty * width + col3] : 0;

		for (int m = 1; m <= (length + 15) / 16; ++m) {
			mArray[ty][tx] = mElement;
			nArray0[ty][tx] = nElement0;
			nArray1[ty][tx] = nElement1;
			nArray2[ty][tx] = nElement2;
			nArray3[ty][tx] = nElement3;
			__syncthreads();

			mElement = (((m * 16 + tx) < length) && (row < height))? matrix1[row * length + m * 16 + tx] : 0;
			nElement0 = (((m * 16 + ty) < length) && (col0 < width))? matrix2[(m * 16 + ty) * width + col0] : 0;
			nElement1 = (((m * 16 + ty) < length) && (col1 < width))? matrix2[(m * 16 + ty) * width + col1] : 0;
			nElement2 = (((m * 16 + ty) < length) && (col2 < width))? matrix2[(m * 16 + ty) * width + col2] : 0;
			nElement3 = (((m * 16 + ty) < length) && (col3 < width))? matrix2[(m * 16 + ty) * width + col3] : 0;
		
			for (int k = 0; k < 8; ++k) {
				value0 += mArray[ty][2 * k] * nArray0[2 * k][tx] +
					mArray[ty][2 * k + 1] * nArray0[2 * k + 1][tx];
				value1 += mArray[ty][2 * k] * nArray1[2 * k][tx] +
					mArray[ty][2 * k + 1] * nArray1[2 * k + 1][tx];
				value2 += mArray[ty][2 * k] * nArray2[2 * k][tx] +
					mArray[ty][2 * k + 1] * nArray2[2 * k + 1][tx];
				value3 += mArray[ty][2 * k] * nArray3[2 * k][tx] +
					mArray[ty][2 * k + 1] * nArray3[2 * k + 1][tx];
			}
			__syncthreads();
		}
		if (row < height) {
			if (col0 < width)
				result[row * width + col0] = value0;
			if (col1 < width)
				result[row * width + col1] = value1;
			if (col2 < width)
				result[row * width + col2] = value2;
			if (col3 < width)
				result[row * width + col3] = value3;
		}
}

__global__  void multiplyKernel_16_prefetch_granularity_1x4_4LoopUnrolling(
	float *result, const float *matrix1, const float *matrix2, const int height, const int width, const int length) { 
		__shared__ float mArray[16][16]; 
		__shared__ float nArray0[16][16];
		__shared__ float nArray1[16][16];
		__shared__ float nArray2[16][16];
		__shared__ float nArray3[16][16];

		int bx = blockIdx.x; 
		int by = blockIdx.y; 
		int tx = threadIdx.x; 
		int ty = threadIdx.y; 

		int row = by * 16 + ty;
		int col0 = 4 * bx * 16 + tx;
		int col1 = (4 * bx + 1) * 16 + tx;
		int col2 = (4 * bx + 2) * 16 + tx;
		int col3 = (4 * bx + 3) * 16 + tx;

		float value0 = 0;
		float value1 = 0;
		float value2 = 0;
		float value3 = 0;

		float mElement = ((tx < length) && (row < height))? matrix1[row * length + tx] : 0;
		float nElement0 = ((ty < length) && (col0 < width))? matrix2[ty * width + col0] : 0;
		float nElement1 = ((ty < length) && (col1 < width))? matrix2[ty * width + col1] : 0;
		float nElement2 = ((ty < length) && (col2 < width))? matrix2[ty * width + col2] : 0;
		float nElement3 = ((ty < length) && (col3 < width))? matrix2[ty * width + col3] : 0;

		for (int m = 1; m <= (length + 15) / 16; ++m) {
			mArray[ty][tx] = mElement;
			nArray0[ty][tx] = nElement0;
			nArray1[ty][tx] = nElement1;
			nArray2[ty][tx] = nElement2;
			nArray3[ty][tx] = nElement3;
			__syncthreads();

			mElement = (((m * 16 + tx) < length) && (row < height))? matrix1[row * length + m * 16 + tx] : 0;
			nElement0 = (((m * 16 + ty) < length) && (col0 < width))? matrix2[(m * 16 + ty) * width + col0] : 0;
			nElement1 = (((m * 16 + ty) < length) && (col1 < width))? matrix2[(m * 16 + ty) * width + col1] : 0;
			nElement2 = (((m * 16 + ty) < length) && (col2 < width))? matrix2[(m * 16 + ty) * width + col2] : 0;
			nElement3 = (((m * 16 + ty) < length) && (col3 < width))? matrix2[(m * 16 + ty) * width + col3] : 0;
		
			for (int k = 0; k < 4; ++k) {
				value0 += mArray[ty][4 * k] * nArray0[4 * k][tx] +
					mArray[ty][4 * k + 1] * nArray0[4 * k + 1][tx] +
					mArray[ty][4 * k + 2] * nArray0[4 * k + 2][tx] +
					mArray[ty][4 * k + 3] * nArray0[4 * k + 3][tx];
				value1 += mArray[ty][4 * k] * nArray1[4 * k][tx] +
					mArray[ty][4 * k + 1] * nArray1[4 * k + 1][tx] +
					mArray[ty][4 * k + 2] * nArray1[4 * k + 2][tx] +
					mArray[ty][4 * k + 3] * nArray1[4 * k + 3][tx];
				value2 += mArray[ty][4 * k] * nArray2[4 * k][tx] +
					mArray[ty][4 * k + 1] * nArray2[4 * k + 1][tx] +
					mArray[ty][4 * k + 2] * nArray2[4 * k + 2][tx] +
					mArray[ty][4 * k + 3] * nArray2[4 * k + 3][tx];
				value3 += mArray[ty][4 * k] * nArray3[4 * k][tx] +
					mArray[ty][4 * k + 1] * nArray3[4 * k + 1][tx] +
					mArray[ty][4 * k + 2] * nArray3[4 * k + 2][tx] +
					mArray[ty][4 * k + 3] * nArray3[4 * k + 3][tx];
			}
			__syncthreads();
		}
		if (row < height) {
			if (col0 < width)
				result[row * width + col0] = value0;
			if (col1 < width)
				result[row * width + col1] = value1;
			if (col2 < width)
				result[row * width + col2] = value2;
			if (col3 < width)
				result[row * width + col3] = value3;
		}
}

__global__  void multiplyKernel_16_prefetch_granularity_1x4_fullLoopUnrolling(
	float *result, const float *matrix1, const float *matrix2, const int height, const int width, const int length) { 
		__shared__ float mArray[16][16]; 
		__shared__ float nArray0[16][16];
		__shared__ float nArray1[16][16];
		__shared__ float nArray2[16][16];
		__shared__ float nArray3[16][16];

		int bx = blockIdx.x; 
		int by = blockIdx.y; 
		int tx = threadIdx.x; 
		int ty = threadIdx.y; 

		int row = by * 16 + ty;
		int col0 = 4 * bx * 16 + tx;
		int col1 = (4 * bx + 1) * 16 + tx;
		int col2 = (4 * bx + 2) * 16 + tx;
		int col3 = (4 * bx + 3) * 16 + tx;

		float value0 = 0;
		float value1 = 0;
		float value2 = 0;
		float value3 = 0;

		float mElement = ((tx < length) && (row < height))? matrix1[row * length + tx] : 0;
		float nElement0 = ((ty < length) && (col0 < width))? matrix2[ty * width + col0] : 0;
		float nElement1 = ((ty < length) && (col1 < width))? matrix2[ty * width + col1] : 0;
		float nElement2 = ((ty < length) && (col2 < width))? matrix2[ty * width + col2] : 0;
		float nElement3 = ((ty < length) && (col3 < width))? matrix2[ty * width + col3] : 0;
		for (int m = 1; m <= (length + 15) / 16; ++m) {
			mArray[ty][tx] = mElement;
			nArray0[ty][tx] = nElement0;
			nArray1[ty][tx] = nElement1;
			nArray2[ty][tx] = nElement2;
			nArray3[ty][tx] = nElement3;
			__syncthreads();

			mElement = (((m * 16 + tx) < length) && (row < height))? matrix1[row * length + m * 16 + tx] : 0;
			nElement0 = (((m * 16 + ty) < length) && (col0 < width))? matrix2[(m * 16 + ty) * width + col0] : 0;
			nElement1 = (((m * 16 + ty) < length) && (col1 < width))? matrix2[(m * 16 + ty) * width + col1] : 0;
			nElement2 = (((m * 16 + ty) < length) && (col2 < width))? matrix2[(m * 16 + ty) * width + col2] : 0;
			nElement3 = (((m * 16 + ty) < length) && (col3 < width))? matrix2[(m * 16 + ty) * width + col3] : 0;
		
			value0 += mArray[ty][0] * nArray0[0][tx] +
				mArray[ty][1] * nArray0[1][tx] +
				mArray[ty][2] * nArray0[2][tx] +
				mArray[ty][3] * nArray0[3][tx] +
				mArray[ty][4] * nArray0[4][tx] +
				mArray[ty][5] * nArray0[5][tx] +
				mArray[ty][6] * nArray0[6][tx] +
				mArray[ty][7] * nArray0[7][tx] +
				mArray[ty][8] * nArray0[8][tx] +
				mArray[ty][9] * nArray0[9][tx] +
				mArray[ty][10] * nArray0[10][tx] +
				mArray[ty][11] * nArray0[11][tx] +
				mArray[ty][12] * nArray0[12][tx] +
				mArray[ty][13] * nArray0[13][tx] +
				mArray[ty][14] * nArray0[14][tx] +
				mArray[ty][15] * nArray0[15][tx];	
			value1 += mArray[ty][0] * nArray1[0][tx] +
				mArray[ty][1] * nArray1[1][tx] +
				mArray[ty][2] * nArray1[2][tx] +
				mArray[ty][3] * nArray1[3][tx] +
				mArray[ty][4] * nArray1[4][tx] +
				mArray[ty][5] * nArray1[5][tx] +
				mArray[ty][6] * nArray1[6][tx] +
				mArray[ty][7] * nArray1[7][tx] +
				mArray[ty][8] * nArray1[8][tx] +
				mArray[ty][9] * nArray1[9][tx] +
				mArray[ty][10] * nArray1[10][tx] +
				mArray[ty][11] * nArray1[11][tx] +
				mArray[ty][12] * nArray1[12][tx] +
				mArray[ty][13] * nArray1[13][tx] +
				mArray[ty][14] * nArray1[14][tx] +
				mArray[ty][15] * nArray1[15][tx];
			value2 += mArray[ty][0] * nArray2[0][tx] +
				mArray[ty][1] * nArray2[1][tx] +
				mArray[ty][2] * nArray2[2][tx] +
				mArray[ty][3] * nArray2[3][tx] +
				mArray[ty][4] * nArray2[4][tx] +
				mArray[ty][5] * nArray2[5][tx] +
				mArray[ty][6] * nArray2[6][tx] +
				mArray[ty][7] * nArray2[7][tx] +
				mArray[ty][8] * nArray2[8][tx] +
				mArray[ty][9] * nArray2[9][tx] +
				mArray[ty][10] * nArray2[10][tx] +
				mArray[ty][11] * nArray2[11][tx] +
				mArray[ty][12] * nArray2[12][tx] +
				mArray[ty][13] * nArray2[13][tx] +
				mArray[ty][14] * nArray2[14][tx] +
				mArray[ty][15] * nArray2[15][tx];
			value3 += mArray[ty][0] * nArray3[0][tx] +
				mArray[ty][1] * nArray3[1][tx] +
				mArray[ty][2] * nArray3[2][tx] +
				mArray[ty][3] * nArray3[3][tx] +
				mArray[ty][4] * nArray3[4][tx] +
				mArray[ty][5] * nArray3[5][tx] +
				mArray[ty][6] * nArray3[6][tx] +
				mArray[ty][7] * nArray3[7][tx] +
				mArray[ty][8] * nArray3[8][tx] +
				mArray[ty][9] * nArray3[9][tx] +
				mArray[ty][10] * nArray3[10][tx] +
				mArray[ty][11] * nArray3[11][tx] +
				mArray[ty][12] * nArray3[12][tx] +
				mArray[ty][13] * nArray3[13][tx] +
				mArray[ty][14] * nArray3[14][tx] +
				mArray[ty][15] * nArray3[15][tx];
			__syncthreads();
		}
		if (row < height) {
			if (col0 < width)
				result[row * width + col0] = value0;
			if (col1 < width)
				result[row * width + col1] = value1;
			if (col2 < width)
				result[row * width + col2] = value2;
			if (col3 < width)
				result[row * width + col3] = value3;
		}
}


